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Abstract 

Recurrent tasks such as pricing, calibration and risk assessment need to be 
executed accurately and in real-time. Simultaneously we observe an increase 
in model sophistication on the one hand and growing demands on the qual¬ 
ity of risk management on the other. To address the resulting computational 
challenges, it is natural to exploit the recurrent nature of these tasks. We 
concentrate on Parametric Option Pricing (POP) and show that polynomial 
interpolation in the parameter space promises to reduce run-times while main¬ 
taining accuracy. The attractive properties of Chebyshev interpolation and its 
tensorized extension enable us to identify criteria for (sub)exponential conver¬ 
gence and explicit error bounds. We show that these results apply to a variety 
of European (basket) options and affine asset models. Numerical experiments 
confirm our findings. Exploring the potential of the method further, we empir¬ 
ically investigate the efficiency of the Chebyshev method for multivariate and 
path-dependent options. 
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1 Introduction 


The development of fast and accurate computational methods for parametric models 
is one of the central issues in computational finance. Financial institutions dedi¬ 
cated to trading or assessment of financial derivatives have to cope with the daily 
tasks of computing numerous characteristic financial quantities. Examples of in¬ 
terest are prices, sensitivities and risk measures for products on different models 
and for varying parameter constellations. With regard to the ever growing mar¬ 
ket activities, more and more of these evaluations need to be delivered in real¬ 
time. In addition we face a const antly rising model so phistication since the original 
work of Black and Scholes ( 197,11 1 and Merton ( 197.'ll h From the early nineties on, 
stochastic volatility and Levy models as well as models based on further classes 
of stochastic processes have been developed that reflect the obse r ved m arket data 
in a more appropriate way. For asset models see e.g. HestonI ( 1993I L Eberlein, 


Keller and Prause (1998), Duffie, Eilipovic and Schachermayer (2003), Cuchiero, 
Keller-Ressel and T e ichma nn (2015). In the case of fixed income models see e.g. 
Eberlein and Ozkan ( 2005l l . Keller-Ressel, Papapantoleon and Teichmann (2013), 


Filipovic, Larsson and Trolle (2014). The aftermath of the financial crisis 2007- 
2009, moreover, has lead to a new generation of more complex models, for instance 
by incorporating more risk factors. The usefulness of a pricing model critically de¬ 
pends on how well it captures the relevant aspects of market reality in its numerical 
implementation. Exploiting new ways to deal with the rising computational com¬ 
plexity therefore supports the evolution of pricing models and touches a core concern 
of present mathematical finance. 

A large body of computational tasks in finance needs to be repeatedly performed 
in real-time for a varying set of parameters. Prominent examples are option pricing 
and hedging of different option sensitivities, e.g. delta and vega, that also need to be 
calculated in real-time. In particular for optimization routines arising in model cal¬ 
ibration, large parameter sets come into play. Purther examples arise in the context 
of risk control and assessment, such as for quantification and monitoring of risk mea¬ 
sures. The following question serves as a starting point of our investigations: How 
to systemically exploit the recurrent nature of parametric computational problems in 
finance with the approached objective to gain efficiency? Looking for answers to this 
question, we focus on Parametric Option Pricing (POP) in the sequel. 

In the present literature on computational methods in finance, complexity reduc¬ 
tion for parametric problems has large ly been addressed by ap plyin g Fouri e r tech - 
niques following the sem inal works of Carr_ancL^Iadan (199^ and R.aibk ( 200(11 1 . 
See also the monograph Boyarchenko and Levendorskiil (2002). Research in this 
area concentrates o n adopting fast Fourier transform (FFT) methods and variants 
for option pricing. Led (2004) accurately describes pricing European options with 
FFT. Further developments are for instance provide d by Lord, Fang, B e rvoet s and 
Poster lee (2008) for early exe r cise o ptions and by Feng and Linetskv ( 2008l l and 
Kudryavtsev and Levendorskii ( 20091 ) for barrier options. Another path to efficiently 


handle large parameter sets that has been pursued in finance relies on reduced basis 
methods. These are te chniques to solve parametrized partial differentia l equat ions. 
Sachs and Schu ( 201oI L Pont, Lantos and Pironneau (2011), Pironneau ( 201 il l and 


Haasdonk, Salomon and Wohlmuth (2012) and Burkovska, Haasdonk, Salomon and 
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Wohlmuth (2015) applied this approach to price European, American plain vanilla 
options and European baskets. EFT methods on the one hand can be highly ben¬ 
eficial when the prices are required in a large number of Fourier variables, e.g. for 
a large set of strikes of European plain vanillas. On the other hand numerical ex¬ 
periments have shown a promising gain in efficiency of reduced basis methods when 
an accurate PDE solver is readily available. In essence all these approaches reveal 
an immense potential of complexity reduction by targeting parameter dependence. 
Hereto, they exploit the functional architecture of the underlying pricing technique 
for varying parameters. 

Financial institutions have to deal simultaneously with a diversity of models, 
a multitude of option types, and—as a consequence—a wide variety of underlying 
pricing techniques. It is therefore worthwhile to explore the possibility of a generic 
complexity reduction method that is independent of the specific pricing technique. 
To do so, we focus on the set of option prices and the set of parameters of interest, 
disregard on purpose the pricing technology and view the option price as a function 
of the parameters. The core idea is now to introduce interpolation of option prices 
in the parameter space as a complexity reduction technique for POP. 

The resulting procedure naturally splits into two phases: Pre-computation and 
real-time evaluation. The first one is also called offline phase while the second is 
also called online phase. In the pre-computation phase the prices are computed 
for some fixed parameter configurations, namely the interpolation nodes. Here, any 
appropriate pricing method—for instance based on Fourier, PDE or even Monte 
Carlo techniques—can be chosen. The real-time evaluation phase then consists of 
the evaluation of the interpolation. Provided that the evaluation of the interpolation 
is faster than the benchmark tool, the scheme permits a gain in efficiency in all cases 
where accuracy can be maintained. Then, we distinguish two use cases: 

• In comparison to the benchmark pricing routine, the fast evaluation of the 
interpolation will eventually outweigh the expensive pre-computation phase, if 
pricing is a task repeatedly employed. Optimization procedures are an obvious 
instance where this feature becomes advantageous. 


• Even if the number of price computations is limited, we can still benefit from 
the separation of the procedure into its two phases. In this way, e.g., idle times 
in the hnancial industry can be put to good use by preparing the interpolation 
for whenever real-time pricing is needed during business activities. 


The question arising at this stage is: Under what circumstances can we hope to find 
an interpolation method that delivers both reliable results and a considerable gain in 
efficiency? 

One could now be tempted to proceed in a naive manner and first define an 
equidistant grid and then interpolate piecewise linearly in the parameter space. 
Numerical experiments for Black&Scholes call prices as function of the volatility, 
for instance, would then provide convincing evidence that the number of nodes 
needed for a given accuracy is considerably high. Increasing the polynomial degree 
might lead to better results. However, convergence might not be guaranteed. iRunge 
(jl90lh showed that polynomial interpolation on equispaced grids may diverge—even 
for analytic functions. Second, the evaluatio n of the poly nomial interpolants may 
be numerically problematic, as it is shown in Rungj ( I90lh that ’’the interpolation 
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problem for polynomial interpolation on a n equidistant gr id is exponentially ill- 
conditioned”, a formulation we borrow from Trefethen ( 201 il l. For these reasons we 
abstain from polynomial interpolation with equidistant grids. Rather we take a step 
back and ask; Which methods for the interpolation of prices as functions of model 
and payoff parameters are numerically promising in terms of convergenee, stability 
and implementational ease? 

Regarding this research question, we need to take into consideration both the 
set of interpolation methods as such and the specific features of the functions we 
investigate. It is well-known that the efficiency of interpolation methods critically 
depends on the degree of regularity of the approximated function. For the core prob¬ 
lem of our study—European (basket) options—we investigate the regularity of the 
option prices as functions of the parameters. We find that these functions are indeed 
analytic for a large set of option types, models and parameters. Taking the perspec¬ 
tive of approximation theory, this inspires the hope to find suitable interpolation 
methods. 

Empirically, we observe that parameters of interest often range within bounded 
intervals. One interpolation method that is highly effective for analytic functions 
on bounded intervals is Chebyshev interpolation. This intensively studied method 
enjoys excellent numerical properties—in stark contrast to polynomial interpolation 
on equally spaced nodal points. The interpolation nodes are known beforehand, im¬ 
plementation is straightforward and the method is numerically stable. Eor univariate 
functions that are several times differentiable, the method converges polynomially 
and for uni variate ana l ytic f unctions convergence is exponential. In a remarkable 
monograph, Trefethen ( 2013l l gives a comprehensive review of Chebyshev interpola¬ 
tion. Rs appealing theoretical properties are indee d of practical use as the soft ware 
tool Chebfur0 demonstrates. In this implemention Platte and Trefethen ( 2008l l aim 
“to combine the feel of symbolics with the speed of numerics”. Therefore Chebyshev 
interpolation is our method of choiceH 

Exploring the potential of interpolation methods for more than one single free pa¬ 
rameter, we choose a tensorized version of Chebyshev interpolation: Eor parameters 
p E where H E IN denotes the dimensionality of the parameter space, the price 
Price^ is approximated by tensorized Chebyshev polynomials Tj with pre-computed 
coefficients cj, j E J, as follows. 


Price^ pa cjTjfp). 


Chebyshev interpolation is a standard numerical method that has proven to 
be extremely useful for applications in such diverse fields as physics, engineering, 
statistics and economics. Nevertheless, for pricing tasks in mathematical finance 
Chebyshe v interpolation still seems to be rarely used and its potential is yet to be 
unfolded. Pistorius and Stolte ( 20121 1 use Chebyshev interpolation of Black&Scholes 
prices in the volatility as an intermediate step to derive a pricing methodology for 


^Chebfun is an open-source software system, see http://www.chebfun.org 

^Chebyshev interpolation shares its good properties wi th for instance Le gendre transformation, 
for which we expect similarly positive results. We refer to iTrefethenI ll201,lh . who states: ”It is the 
clustering near the ends of the interval that makes the difference, and other sets of points with 
similar clustering, like Legendre points [...] have similarly good behaviour.” 


4 


































a time-changed model. Independently from us, Pachon ( 20161 1 recently proposed 
Chebyshev interpolation as a quadrature rule for the computation of option prices 
with a Fourier type representation, which is comparable to the cosine method. 

Our main results are the following: 


• Theorem 13.21 provides accessible sufficient conditions on options and models 
that guarantee an asymptotic error decay of order 0(^g~ in the total 
number N of interpolation nodes where > 1 is given by the domain of 
analyticity and D is the number of varying parameters. 


• More specific conditions for parametric European options in Levy models are 
provided in Corollarv l3.6[ while Corollary 13.91 provides the framework for para¬ 
metric basket options in affine models. 


These results establish an error analysis that is based on the domain of analyticity 
of the prices as functions of the parameters. Observing that typical payoff functions 
are not smooth, we cannot expect an exponential error decay for interpolation with 
arbitrarily small maturities. Small maturities thus serve as an example that domains 
of analyticity need to be carefully studied. 

• The investigations in Sections I3.2H4.2I show that for a large set of relevant 
(basket) options, models and free parameters a domain of analyticity can in¬ 
deed be identified. This gives examples of relevant financial applications where 
(sub)exponential error decay is guaranteed. 

To numerically validate the theoretical results we compare prices obtained by Cheby¬ 
shev interpolation to benchmark prices by Fourier techniques. 


• Numerical experiments in affine models confirm the theoretical error decay em¬ 
pirically for European call options (Figure 15.3p and digital down&out options 
(Figure [531) • the considered model examples of Black&Scholes, Merton, 
CGMY and Heston we observe L°°-error levels of order 10“^° using not more 
than N = (25 -|- 1)^ Chebyshev interpolation nodes when D = 2 parameters 
are varied. 


Numerical results show that already a small number of nodes leads to high accuracy. 
This motivates us to further explore the potential of the Chebyshev method for 
multivariate options. Here we deliberately go beyond the scope of our theoretical 
results and consider additional features like path-dependency. 

• For multivariate basket and path-dependent options in the Black&Scholes, 
Heston and Merton model we use Monte Carlo as reference method. In all 
of our settings in Section 15.21 Chebyshev interpolation achieves an accuracy 
that is similar to the accuracy of the Monte Carlo simulation itself (10“^) for 
D = 2. 


In addition we present empirical results demonstrating the efficiency of the Cheby¬ 
shev method. 

• The gain in efficiency in comparison to Fourier techniques is first validated for 
bivariate options of European type in Section 15.3.11 
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• Secondly, the explicit gains in efficiency in comparison to Monte Carlo methods 
are shown in Section 15.3.21 taking multivariate lookback options in the Heston 
model as examples. 

The remainder of the article is organized as follows. In Section [2] we intro¬ 
duce Chebyshev interpolation in detail and present the general convergence results. 
Section [3] establishes a convergence analysis of Chebyshev interpolation for POP. 
We formulate analyticity conditions for the payoff profiles and models that guaran¬ 
tee (sub)exponential convergence of the method. These conditions are verified in 
Section [3] for different option types, models and free parameters. The numerical 
experiments in Section 5 confirm these findings using Fourier techniques. Pricing 
basket options, the gain in efficiency is numerically investigated. Experiments based 
on Monte Carlo and finite differences moreover suggest to further explore the po¬ 
tential of the approach beyond the scope of the theoretical investigations from the 
previous sections. The resulting conclusion and outlook are presented in Section 6. 
Finally, the appendix provides the proof of the multivariate convergence result. 


2 Chebyshev Polynomial Interpolation 


In this se c tion w e introduce the notation for Chebyshev interpolation. Following 
Trefethen ( 20131 1. the one-dimensional version is shown. Then we present the mul¬ 


tivariate extension and convergence results. Consider an option price with a single 
varying parameter 


(2.1) Price^, pG[—1,1]. 

An interpolation of Price^ with Chebyshev polynomials of degree N is of the form 

N 

(2.2) I]\f{Price^'^){p) := 

j=0 


with coefficients 


(2.3) 


c, := cos 

k=0 


N 


j<N, 


and basis functions 


(2.4) Tj{p) := cos (j arccos(p)) for p G [—1,1] and j < N 


where ^ indicates that the first and last summands are halved. The Chebyshev 
nodes pk = cos (vr^) may conveniently be displayed in a graph, see Figure 12.11 
For an arbitrary compact parameter interval [p,p], interpolation (12.211 needs to be 
adjusted by the appropriate linear transformation. 
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Chebyshev nodes for D=1 



Figure 2.1: A set of Chebyshev points Pk € [—1,1] (blue) for degree = 20 and equidistantly 
spaced auxiliary construction points (red) on the semi-circle. 


2.1 Multivariate Chebyshev Interpolation 


The Chebyshev polynomial in terpolation (I2.2D-(I2.4) has a tensor based extension to 
the multivariate case, see e.g. Sauter and Schwab ( 20041 ) . In order to obtain a nice 
notation, consider the interpolation of prices 


(2.5) 


Price^, pG[—1,1]^. 


For a more general hyperrectangular parameter space V = \p^,p^]x .. .x \p^,p^], the 
appropriate linear transformations need to be performed. Let N := (A^i,... ,N£)) 
with Ni G INq for i = 1,... ,D. The interpolation with + 1) summands is 

given by 

(2.6) %(Phce('))(p) :='^CjTj{p), 

3&J 

where the summation index j is a multiindex ranging over J := {(ji,... ^jo) £ ■ 

ji < A^ifor i = 1,... ,D}, i.e. 

N\ No 

(2.7) /^(Phce(-))(p) =Y1 Y. C(iiv..,m)%,...,m)W- 

ji=o iD=o 

The basis functions Tj for j = {ji, ■ ■ ■, jo) £ J are defined by 

D 

(2.8) Tj{pi,... ,pd) = '[lTj.{pi). 

i=l 


The coefficients Cj for j = {ji, ■ ■ ■, jo) £ J are given by 


(2.9) 


(n 




2=1 


Ni 


No D 

{ki,...,ko) 


^"... "PnceP' 

fci=0 k]j=0 


n 

2=1 


COS JiTT 


, ^ ki 
Ni 
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where ^ " indicates that the first and last summand are halved and the Chebyshev 
nodes for multiindex k = (fci ,... G J are given by 

(2.10) p’" = {Pki,-■ ■ ,Pkn) 

with the univariate Chebyshev nodes = cos ~ 0,... ,iVj and i = 

1,... ,D. A set of D-variate Chebyshev nodes for D = 2 and A^i = A ^2 = 20 

is depicted in Figure [22J 
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Chebyshev nodes for D=2 



Ph 


Figure 2.2: A set of D-variate Chebyshev points p* £ [—1,1]^ for D = 2 and A^i = N 2 = 20. 


2.2 Convergence of Multivariate Chebyshev Interpolation 

In the univariate case, it is well known that the error of approximation with Cheby¬ 
shev polynomials decays polynomially for differentiable functions and exponentially 
for analytic functions. Let / be analytic in [—1,1] then it has an analytic extension 
to some Bernstein ellipse i?([—l,l],p) with parameter p > 1, dehned as the open 
region in the complex plane bounded by the ellipse with foci ±1 and semiminor and 
semimajor axis lengths summing up to q. This and the following result traces back 
to the seminal work of Bernstein ( 1912I L 


Theorem 2.1. 


Trefethejl \ 20 ni. Theorem 8.2) Let a function f be analytic in the 
open Bernstein ellipse B([—1,1], p), with g > 1, where it satisfies \ f\ < V for some 
constant C > 0. Then for each N > 0, 


\\f-lN{f)\\L^il-l,l])<^V 


..-N 


p - 1 
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In the multivariate case we will extend a convergence result from 
(|2004l i. We consider parametric option prices of form 


Sauter and Schwab 


(2.11) Price^ for p G P 

with V C IR^ of hyperrectangular structure, i.e. V = \p^,p^]x .. .x [p^,p^] with real 
p. < p. for all i = 1,... ,D. We define the D-variate and transformed analogon of 
a Bernstein ellipse around the hyperrectangle V with parameter vector g € (l,oo)^ 
as 


(2.12) B{V, g) := B{\p.^,p^], gi) x ... x B{\p^,p^], go) 

with B{\p,]i\,g) := Tjppj o i?([—1,1],^), where for p G C we have the transform 

%?)] (3^(p)) := P + ^(l - ^{p)) and r[p^p](9(p)) := ^’^{p). We call B{'P,g) 
generalized Bernstein ellipse if the sets B{[—1,1], gi) are Bernstein ellipses for i = 

1,...,Z1. 

Theorem 2.2. Let V 3 p Priced he a real valued function that has an analytic 
extension to some generalized Bernstein ellipse B{'P,g) for some parameter vector 
g G (l,oo)^ and sup^g^^-p |Price^| < V. Then 

! D 

max |Phce^ — /■^(Phce^'^)(p)| < B • j ^ -z 

\i=i j=i ^ ~ 

The proof of the theorem is provided in Appendix 1X1 

Corollary 2.3. Under the assumptions of Theorem \2.S\ there exists a constant C > 0 
such that 



(2.13) 


max 

P&V 


\Price^ — Ij^{Price^'^){p)\ < Cg —, 


where g = min u,- and N = min W. 

- i<i<D^ — l<i<D 

Remark 2.4. In particular, under the assumptions required by Theorem \2.2\ with 
N = n?=l{N^ + 1) denoting the total number of nodes, Corollary \ 2. .shows that the 
error decay is of (sub)exponential order 0[g~ for some g > 1. 


In the setting of Theorem 12.21 additionally the derivatives of Price^ are ap¬ 
proximated by the according derivativ e s of t he Chebyshev interpolation. The one- 
dimensional result is shown in Tadmor ( IhSfil i and a multivariate result is derived i 


m 


Canuto and Quarteroni (1982) for functions in Sobolev spaces. These results allow 


us to obtain the Chebyshev approximation of deriv atives with no additional cost . 
To state the according convergence results we follow ICanuto and Quarteronil (jl982l l 
and introduce the weighted Sobolev spaces for u G M by 


(2.14) = {cfG L\V) : UWwnv) < oo} , 


with norm 





















wherein a = (ai, ...,a£))€M^isa multiindex and 5 “ = 9 “^ • • • ( 9 "® and the weight 
function ui on V given by 


D 


ivix) := H j(x,-)), cu(r[/ ](x,-)) := (1 -^ 

—3 —3 —3 


i=i 


_ p■~p / \ 

with Tjp p,] (p) = + . Then we apply the result of iCanuto and Quarteroni 

(|1982l . Theorem 3.1) in the following corollary. 

Corollary 2.5. Under the assumptions of Theorem \2.^ for any y < cr E IN and 
any cr > /r E INq there exists a eonstant C > 0 such that 

||Pnce(') < CN‘^^-^\\Price^^^\\^^,^^Jyy 

Proof. In our setting we have V C of hyperrectangular structure. Under the 
assumptions of Theorem 12.21 it follows th at p i— Price^ E Wf(V ) and therewith 
p I—)• Price^ E Wf’‘^{V). Before we apply Canute and Quarteronil ( 19821 . Theorem 
3.1), which assumes V = [—1,1]'^, we investigate how the linear transformation r-p, 
as introduced in the proof of Theorem l2.2l influences the derivatives. Let p i—)• Price^ 
be a function on V. We set h{p) = Price^ o Tp{p). Further, let Ij^{h){p) be the 
Chebyshev interpolation of h{p) on [—1,1]^. Then, it directly follows 

PriceP - %(Phce('))(p) = (h{-) - T-^(h){-)^ o T^^(p). 

First, let us assume D = 1, i.e. V = \p,p], and let a E Nq. For the partial derivatives 
it holds 


d°‘PriceP — d°‘Ij^{Price^'^){p) = 9 “ yPriceP — Ij^{Price^'^ 

= ((^(O - %(^)(-)) o 

= {d^h{Tf;\p)) - d%{h^-)){Tf;\p))) 


= & 


A-l 


P — P 


d^h {t^\p))- (rp^(p))). 


Repeating this step iteratively yields 
d°‘PriceP — d°‘IjT(Price^'^){p) = - 


(p - pY 


d^h 


(TPip)) - {tP{p[ 


Therewith, the error on [— 1 , 1 ] is scaled with a factor Extending this to 

the D-variate case with, this analogously results with a = (ai,..., an) E IN^ is a 
multiindex and 9 “ = 9 “^ • • • 9 “^ 


d'^PriceP — 8°^ Ij^{Price^Y{p) — 

^ 2\°‘i\ 


n 


{Pi - Pi)I"'I 


9 “h 


{TpHp))- 9“Iiv(h(')) (t^Hp))) 
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From Theorem 3.1 in 
for h{-) oillV= [—1,1]^, i.e. for any y < cr G IN and any a > fj, G Mq there exists a 
constant C > 0 such that 

(2.16) ||h(-) - (1^)’ 

For arbitrary V the constant from l2.16l has to be multiplied with the according factor 
resulting from the linear transformation r-p. □ 


Canuto and Quarteronil ( 198211 the assertion follows directly 


The result in Corollary 12.51 is given in terms of weighted Sobolev norms. In the 
following remark we connect the approximation error in the weighted Sobolev norm 
to the C\V) norm, where C\V) is the Banach space of all functions u in V such 
that u and d°^u with |a| < / are uniformly continuous in V and the norm 


is finite. 


= maxmax |5 “m(p)| 
^ ’ |a|<Z P&V 


Corollary 2.6. Under the assumptions of Theorem \2.i^ for any y < cr G IN and any 
a > pL G INo and I G INq with n — I > ^, there exists a constant C{a) > 0 depending 
on a, such that 

||Przce^'^ — Ijq{Price max sup 19“ Price^ \. 

l"l<o'pe'P 


Proof. In the setting of Corollary 12.51 we start with the estimation of the approxi¬ 
mation error in the weighted Sobolev norms, 

(2.17) ||Pnce» - %(Phce(-))(-)||u^.,^(p) < CAr2/^-"||Pnce« 

On V it holds that w{p) > 1 and therewith we can deduce for the Sobolev norm 
with a constant weight of 1, 

||Price(') -/^(Phce('))(-)||u/2^.“(p) > UPrice^') -/^(Price('))(-)||^M,i(p). 

With W^iV) the usual Sobolev space. 


= {<(>£ : U\\wf{r) < oo} ’ M\wf(V) = X] / 


Corollary 6.2 from Wloka ( 1987l i directly yields that for any I with fj, — I > ^ there 
exists a constant C such that 


(2.18) ||Phce(') - %(Pnce('))(-)||ch7^) ^ CWPrice'^'^ - Ijp{Price^'^){-)\\w^(jyy 

In formula (I2.18P we have derived a lower bound for the left hand side of expression 
(I2.17p . Next, we will find an upper bound for the right hand side of (|2.17p . Prom 
the definition of the weighted Sobolev norm, see (|2.15p . it follows 


II Price I 




I |(9"Price^pa;(p) dp 


; <cr 


Ip 


< 


yy sup |5"Price^p / a;(p) dp. 


\ \a\<o-P^P 


P 
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Here, we apply fpUj{p)dp = and that there exists a constant a'^(a) depending 
on a such that 



(2.19) 


= a{a) max sup |9“Phce^|7r 2 . 
l«l<o'pe'P 


Finally, using ()2.18p and (I2.19p in (I2.17p yields an estimate of the approximation 


error in the C\V) norm, 



Collecting all constants in C{a), we achieve 

llPrice^'^ — max sup |9"Phce^|. 


3 Exponential Convergence of Chebyshev 
Interpolation for POP 

In this section we embed the multivariate Chebyshev interpolation into the option 
pricing framework. We provide sufficient conditions under which option prices de¬ 
pend analytically on the parameters. 

Analytic properties of option prices can be conveniently studied in terms of 
Fourier transforms. First, Fourier representations of option prices are explicitly 
available for a large class of both option types and asset models. Second, Fourier 
transformation unveils the analytic properties of both the payoff structure and the 
distribution of the underlying stochastic quantity in a beautiful way. By contrast, if 
option prices are represented as expectations, their analyticity in the parameters is 
hidden. For example the function K 1 —)• (St — K)^ is not even differentiable, whereas 
the Fourier transform of the dampened call payoff function evidently is analytic in 
the strike, compare Table 14.11 on page [T9l 

3.1 Conditions for Exponential Convergence 

Let us first introduce a general option pricing framework. We consider option prices 
of the form 


(3.1) 


PriceP=^P"’P"^ = E{fP\xP^)) 


where fP^ is a parametrized family of measurable payoff functions fP^ : —)• 1R+ 
with payoff parameters p^ and XP is a family of IR'^-valued random variables 
with model parameters p'^ € The parameter set 


(3.2) 


p = {p\p‘^) eV = V^ xV^ 


12 









is again of hyperrectangular structure, i.e. x ... x and = 

\Pra+v^m+^ X ... X for some 1 < m < D and real < p. for all i = 1,..., D. 

Typically we are given a parametrized 11“^-valued driving stochastic process 
with being the vector of asset price processes modeled as an exponential of , 

i.e. 

(3.3) Sf = Sq exp(i/f ’*), 0 <t <T, 1 < i < d, 

and is an J^y-measurable R'^-valued random variable, possibly depending on 
the history of the d driving processes, i.e. p^ = (T,p') and 

XP^ := 0 < t < T), 

where ih is an R'^-valued measurable functional. 

We now focus on the case that the price (I3T]) is given in terms of Fourier trans¬ 
forms. This enables us to provide sufficient conditions under which the parametrized 
prices have an analytic extension to an appropriate generalized Bernstein ellipse. For 
most relevant options, the payoff profile fP is not integrable and its Fourier trans¬ 
form over the real axis is not well defined. Instead, there exists an exponential 
dampening factor rj G R*^ such that fP G L^(R'^). We therefore introduce 
exponential weights in our set of conditions and denote the Fourier transform of 
g G L^(R‘^) by 


g{z) := [ g{x) dx 

Jud 

and we denote the Fourier transform of / G L^(R'^) by f{-—ig). The exponential 

weight of the payoff will be compensated by exponentially weighting the distribution 
2 2 
of XP and that weight will reappear in the argument of pP , the characteristic 

function of XP^ . 

Conditions 3.1. Let parameter set V = x R^ of hyperrectangular 

structure as in (I 321 ). Let Q G (l,oo)^ and denote := {Qi,---,gm) and g^ := 
{Qm+i, ■ ■ ■, Qd) and let weight g G R*^. 

(Al) For every p^ G the mapping x fP^{x) is in L^(R'^). 

(A2) For every z G R'^ the mapping p^ i— fP'^iz — ig) is analytic in the general¬ 
ized Bernstein ellipse B{V^, g^) and there are constants ci,C 2 > 0 such that 
suppigB(-pi^gi) \fP^{—z — ig)\ < cie'^^bl for all 2: G R'^. 

2 

(A3) For every p^ G the exponential moment condition E(^e~^P’^^ ^) < 00 
holds. 

(A4) For every 2; G R'^ the mapping 1—>■ ipP^ {z -|- ig) is analytic in the generalized 

Bernstein ellipse B{V‘^,g‘^) and there are constants a G (1,2] and ci,C2 > 0 
such that supp 2 gg(-p 2 ^g 2 ) \(pP {z + ig)\ < ci for all 2; G R'^. 

Conditions (Al)-(A4) are satisfied for a large class of payoff functions and asset 
models, see Sections O and 01 
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Theorem 3.2. Let g G (l,oo)^ and weight rj € Under conditions (Al)-(A4) 
we have 


max 

p£V 


\Price^ — Ij^{Price^'^){p)\ 


D 

< 2^+^ • V ■ 


D 


D 




i=l 


i=i 


1- g 


-2 


1 

2 


The proof of the theorem builds on the following proposition that can be derived 
from Eberlein, Gian and Papapantoleon (2010, Theorem 3.2), observing that the 

proof solely uses that z fP^{—z — iri)ipP‘^ {z + iri) belongs to L^(R'^) instead of 
the slightly stronger assumption that z gP (-^ + iv) belongs to L^(R'^) which is 
imposed there. 

Proposition 3.3. Assume conditions (Al), (A3) and that the mapping z 
ir])(pP (z + irf) belongs to L^(R'^) for every p = € V, then the option 

prices (1311 have for every p = {p^,p^) G V the Fourier transform based repre¬ 
sentation 

(3.4) PriceP =[ fP\-z)ipP\z) dz. 

(^Tt) 

We are now in a position to prove Theorem 13.21 

Proof of Theorem 1,9.111 In view of Theorem 12.21 it suffices to show that the mapping 
p I—)• PriceP has an analytic extension to the generalized Bernstein ellipse B[V,g). 
Thanks to Proposition [TSl we have the following Fourier representation of the option 
prices. 


= 7T~fd f fP\-z)<pP\z) dz. 

(^TT) 

Due to assumptions (A2) and (A4) the mapping 

P = (p\p^) P\-z)tp\z) 


has an analytic extension to B{V, g). Let 7 be a contour of a compact triangle in the 
interior of B{\p.,p.], gi) for arbitrary i = 1,..., D. Then thanks to assumption (A2) 
and (A4) we may apply Fubini’s theorem and obtain 


Price^P^''"'P°\z) dpi = 


1 


(27r) J -f J 'U3‘-\-iri 

1 r f 


(271") JIRd+ir; . 


fp\-z)(pP^ {z) dzdpi 
fP^ {-z)ipP^ (z) dpi dz = 0 . 


Moreover, thanks to assumptions (A2) and (A4), dominated convergence shows 
continuity of p 1 —>■ PriceP in B(V,g) which yields the analytic ity of P 7 — )• P riceP 
in B{V, g) thanks to a version of Morera’s theorem provided in iJanichl (j2004l . Satz 
8 ). □ 
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Similar to Corollary 12.51 in the setting of Theorem l3.2l additionally the according 
deriyatiyes are approximated as well by the Chebyshey interpolation. A yery in¬ 
teresting application of this result in finance is the computation of sensitiyities like 
delta or vega of an option price for risk assessment purposes. Theorem 13.21 together 
with Corollary 12.51 yield the following corollary. 

Corollary 3.4. Under the assumptions of Theorem \3.‘A for all Z G N, /r and a with 
a>^,0<p<a and p — I > ^ there exist a constant C, such that 

\\PriceP - Ij^{Price’^'"i){p)\\ci(j,) < CNff^^~^\\PriceP\\w^(j,), 
where the spaces and norms are defined in Section fKB . 


Remark 3.5. The upper bound in condition (A2) is tailored to the application to 
payoff profiles with varying option parameters, compare Section \4fl\ and particularly 
Table \4fl\ below. Examples of (time-inhomogeneous) Levy yro cesses satisf ying the 


upper bound in (Af) for = [T,T] C ( 0 , oo) are provided in GlaA i '201f\) . where 


also its connection to the parabolicity of the Kolmogorov equation of t he process i s 
investigated. In a fix parameter setting, condition (A3) for a G (1,2] in GIoaI \201^) 
translates to our upper bound in (A4). This condition is satisfied for a large variety 
of models. 


3.2 Convergence Results for Selected Option Prices 

In the preyious subsection, Conditions 13.II and Theorem 13.21 introduced a framework 
in which the Chebyshey approximation achieyes (sub)exponential error decay. Now 
we relate this abstract framework to two specific settings. 


3.2.1 European Options in Univariate Levy Models 

Let r be the deterministic and constant interest rate. We consider the parametrized 
family of asset prices, 

(3.5) Sf := So 


with t > 0. For fixed vr = (a, b) let be a Leyy process with characteristics 
{a‘^,b,F) and |x|F(dx) < oo. Due to the Theorem of Leyy-Khintchine we 

have the following representation of the characteristic function of the parametrized 
Levy process 

(3.6) E(e*^'^") =-^ + ibz + [ -1 - izx)F{dx). 

2 Jr 

We assume is defined under a risk neutral measure. Therefore, for every vr G ft 
we assume ) < oo for some and equivalently all t > 0 and the drift condition 

(3.7) b = b{r,a) = r — ^— / (e* —1 — x)F(dx), 

2 Jr 

to ensure that the discounted asset price process is a martingale. Denoting 


(3.8) 


filz) := f —1 — izx)F{dx), 

Jr 
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the pure jump part of the exponent of the characteristic function, this reads as 

2 

(3.9) 5 = 6 (r,o-) = r-y 

The fair value at time t = 0 of a European option with payoff function for 
K := [I£, iiT] C R with maturity T € [T, T] C (0, oo) is given by 

(3.10) Price = e"’’^ P(/^(So e^^)). 

In order to guarantee (sub)exponential convergence of the Chebyshev interpolation, 
in the following corollary, we translate the exponential moment condition (A3), the 
analyticity condition and the upper bound in (A4) to conditions on the cumulant 
generating function . 

Corollary 3.6. Let Conditions (Al) and \(A2^ he satisfied for weight r/ G R and 
0 ^ > 1 and set = [A, A]. Moreover, let = [S'oj'S'o] x [T_,T] C R^ with 

•So) E > 0. Assume 


(3.11) 




(e-^x ve^)F(dx) < oo 


and there exist constants Ci, 6*2 >0 and fi <2 such that 
(3.12) ir])) <Ci + C 2 \z\^ for all z 

If additionally one of the following conditions is satisfied, 

(i) u > 0 , 


(a) there exist a G (1,2] with fi < a and constants Ci,C 2 > 0 such that 
iR-[fi){z + irj) < Cl — C2\z\‘^ forallzGM., 


then there exist constants C > 0 and g > 1 such that 

(3.13) max ^ \ Price^ — I-^{Price^'^){p)\ < Cq~—, 


where N = min A 3 . 

l<i <3 


Proof. In view of Theorem 13.21 and Corollary 12.31 under the hypothesis of this 
corollary, it suffices to verify Con ditions (A3 ) and (A4). Thanks to the expo¬ 


nential moment condition (|3.11l) . Sato ( 19991 . Theorem 25.17) implies that the 


validity of the Levy-Khintchine formula (13.6p extends to the complex strip U := 
R -|- i ([0,1] U (sgn(r/)[0, |??|])), i.e. for every tt € 11 and every t > 0, 


(3.14) 


E (e“^" ) = for every zGU. 


In particular, A(e < 00 for every vr G H, which yields (A3). Denote 

So := log(So)- The assumptions of the corollary yield for every z G U the analyticity 


of 

(so,t) ^ 
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on the generalized Bernstein ellipse g^) for some parameter G (1 ,oo)^, i.e. 

the validity of the analyticity Condition in (A4). Next we consider 

\^'P\z + iv)\<\ e“o(^+*^) 11 I 

< exp {—SQf? — 9‘(so)^ + {z + ir])) — {z + irj ))} . 

An elementary manipulation shows for every z G R, 

cr^z'^ r 

{z + it])) ='ij;^{if]) ---h / (cos(zx) —l)e ^*F(dx) 

2 Jr 

which in combination with the assumption on the imaginary part of ^|J yields the 
upper bound in Condition (A4) under assumption (ii). In view of 


IR 


(cos( 2 ;x) — l) e F{dx) < 0, 


condition (A4) also follows under assumption (i). 


□ 


Let us notice that for the Merton model, Condition (|3.11l) is satisfied for ev¬ 
ery r/ G R and Condition (I3.12p is satisfied with /3 = 1. Examples of Levy jump 
diffusion modes, i.e. examples satisfying condition (i) of Corollary 13.61 are e.g. the 
Black&Scholes and the Merton model. Examples of p ure jump Le vy models sat¬ 
isfying condition (ii) of Corollary 13.61 are provided in Glau ( 20161 1. compare also 
Remark 13.51 and Table 141^ below. 

Eor simplicity and in accordance with our numerical experiments, we fixed the 
model parameters tt. The analysis can be extended to varying model parameters vr. 


Remark 3.7. Assuming (|3.11l) . we observe that the mappings 


2 ^ 

(3.15) (r,( t) 1 -^ 6(r,( t) = r — — / (e^ —1 — x)E(dx) 

2 Jr 

and {r,a) i-G (z) for every z ^ U are holomorphie. Moreover, {r,a) e-)• 

b{r,a) is bounded on every generalized Bernstein ellipse B{[r,r] x [a,a],p) with 
g G (l,oo)2. 

Remark 3.8. Corollary 13.61 allows us to determine explicit error bounds for call 
options in Levy models. The fair price at t = 0 of a call option with strike K and 
maturity T with deterministic interest rate r >0 is 

(3.16) Calll°'^ = E{So -K)^ 


under a risk-neutral probability measure. Noticing that 

(3.17) Call^°’^ = KE{{So/K)e^^ -l) + , 


it suffices to interpolate the function 

(3.18) {T,So)r^ Cal4°’^ 
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on [T_,T] X \So / K, Sp/K ] in order to approximate the prices CaU^°'^ for values 
{T,So,K) G [T, T] X [5*0 , 5o] x [^,K] C (0,oo)^. This effectively reduces the dimen¬ 
sionality D of the interpolation problem by one. 

Let us fix some rj < —I and let V = \Ti,T] x \Sq/ K, Sq/ K ], o,nd 

Cf = =^, then for every G (1, C,^ + \/(C^T “ 1) for j = 1, 2, let q = {g^, g^) and 


V := sup 

(T,So)GB(P,^) 


Cal4°'^ , 


then 


max ICall^’^ 
(T,So)eP ' 


lN,MCall\]’^){T,So)\<W 




In particular, there exists a constant C > 0 such that 

(3.19) max \Cald°’^ - In, No{Cal6'V)(T, So)\ < Cg-^, 


where g = mm{^)i, g 2 } and N_ = minjiVi, A^ 2 }- 
Additionally, when fixing the maturity T, letting ( 
nential error decay 


SoK+^K 

sEk-^’ 


we obtain the expo- 


(3.20) 


_max_ I Call^'^ 

^/K<So<So/K 


I^{Call^^^’^){So)\<4vj^ 


for some g G 


(1,C +V? 


1) and V 


sup 

SoeB([^/K;sE/K],g) 


Call^°'^ 


3.2.2 Basket Options in AfRne Models 

Let X'"' be a parametric family of affine processes with state space D C IR,'^ for 
tt' G n' such that for every vr' G 11' there exists a complex-valued function and a 
C'^-valued function cjT such that 


(3.21) 


T 


P^ = {t,xy) {t,iz)+{y' {t,iz),x) ^ 


for every t > 0, z ^ and x G S. Under mild regularity conditions, the functions 
y and (jT are determined as solutions to generalized Riccati equations. We refer 
to DufRe, Filipovic and Schachermayer (2003) for a detailed exposition. The rich 
class of affine processes comprises the class of Levy processes, for which y (t, iz) = 
tif^ (z) with given as some exponent in the Levy-Khintchine formula (13.6p and 
(jT [t, iz) = 0. Moreover, many popular stochastic volatility models such as the 
H eston model as well as stochastic vola tility models with jumps, e.g. the model 
of iBarndorff-Nielsen and Shephardl (1200 lli and time -changed Levy models, see Carr, 
Geman, Madan and Yor (2003) and Kail sen ( 2006l i. are driven by affine processes. 

Consider option prices of the form 


(3.22) 


Price = E{f^{Xf')\XS' = x) 


where f^ is a parametrized family of measurable payoff functions f^ : 
for K G Vf 


IR-i 
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Corollary 3.9. Under the conditions (Al)-(A3) for weight rj € g G (l,oo)'^ 
and V = V^ <Z of hyperrectangular structure assume 

(i) for every parameter = (t,x, 7 r') gV'^C R^“™ that the validity of the affine 
property (| 3 . 2 ip extends to z = Ji + irj, i.e. for every z G R + ir/, 

|xj' = x) = ihiz)+{r'it,iz),x)^ 

(a) for every z G R'^ that the mappings (t, vr') i—)• {t,iz — p) and {t,Tr') 

(jf (t,iz — p) have an analytic extension to the Bernstein ellipse BfU', g') for 
some parameter g' G (1,00)^“™“^, 

(Hi) there exist a G ( 1 , 2 ] and constants Ci,C2 > 0 such that uniformly in the 
parameters p^ = {t,x,Tr') G B{V^,ffi) for some generalized Bernstein ellipse 
with ^ G (1,00)^“™ 

3 ft(t, iz — rf) + {(ffi {t, iz — r]),x)) < Ci — C'2|2;|“ for all 2; G R. 


Then there exist constants C > 0, g > 1 such that 


max ^ |Price^ — Ij^{Price^'^){p)\ < Cg 

Proof Thanks to Theorem 13.21 and Corollary 12.31 and in view of the assumed validity 
of Conditions (A1)~(A3), it suffices to verify Condition (A4). While assumptions 
(i) and (ii) together yield the analyticity condition in (A4), part (hi) provides the 
upper bound in (A4). □ 


The existence of exponenti a l mom ents of affine processes has been investigated by 
Keller-Ressel and Maverhofer ( 20151 1. where also criteria are provided under which 


formula (I3.2ip and the related generalized Riccati system can be extended to complex 
exponential moments z G The question has already b een treated for important 
specia l cases, which allow for mor e explicit conditions. [Filipovic and Maverhofer 


( 20091 1 consider affine diffusions and Cheridito and Wugalter ( 2012 1 investigate affine 
processes with killing when the jump measures possess exponential moments of all 
orders. 


4 Analyticity Properties 

4.1 Analyticity Properties of Selected Payoff Profiles 

We first list in Table 14.11 a selection of payoff profiles f^ for option parameter K 
as function of the logarithm of the underlying. As we have seen in Proposition 13.31 
we can represent option prices under certain conditions by their Fourier transform. 
Therefore, the table provides the Fourier transform f^ of the respective option 
payoff, as well. All examples in Table 14.11 are special cases by k := log(A) of the 
following lemma that has a straightforward proof. 
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Type 

Payoff 

Weight 

Fourier transform 

f^ holomor- 


fix) 

V 

1 

. 

phic in log(i^) 

Call 

(e^ -K)+ 

< -1 

j^iz + l + r] 

{iz+ri){iz+l+ri) 

/ 

Put 

{K - e^)+ 

> 0 

j^iz + l + r} 

{iz+ri){iz+l+ri) 

/ 

Digital 

down&out 

^a:>log(i^) 

< 0 

J^iz + 7] 

iz-\-r} 

/ 

Asset-or-nothing 
down&out 


< -1 

^iz-\-l-\-rj 

/ 


Table 4.1: Examples of payout profiles of a single underlying. 


Proposition 4.1. Let K, x 9 {k^x) i-T- f^{x) he of form f^{x) = h{k) f^ {r^{x)) 
with a transform T^{x)j = Xj + Ujk with aj G R for every j = and a 

holomorphic function h. Let ij G such that x f^{x) belongs to L^(R'^) 

for every A: G R. Then 

(i) the mapping k f^{z — irj) has a holomorphic extension and 

(ii) \f^{z — iri)\ < \h{k)\ for every fc G C. 

Proof. Both assertions are immediate consequences of 

f{T^){z-ir]) =expi -k'^aj{izj - rj) >f{z-id) 

^ i=i ^ 

for all G R*^ and all k G C. □ 

Example 4.2 (Call on the minimum of several assets). The payoff profile of a call 
option on the minimum of d assets is defined as 

(4.1) f^{x) = (e^i A A • • • A - K)+ , 

for X = {xi,... ,Xd) G R'^ and strike K G R'*". With dampening constant rj G 
(— oo, — 1 )“^ we have x eA f^{x) G L^(R'^). Proposition \4-l\ shows that K i-A 
f^{z — irj) has a holomorphic extension. 

For some payoff profiles it is not possible to transform them to integrable func¬ 
tions by a multiplication with an exponential dampening factor. Thus we split them 
into summands which can be appropriately dampened. Therefore, we decompose 
the unity, 1 = Y^j=i (^) with the distinct orthants of R'^. More precisely, 
for j = 1,..., 2“^ let := (Ci, • • •, C^) with c/ G {—1,1} for the 2‘^ different possible 
configurations and let 

(4.2) := {(xi,..., xf) G R"^ | (jXi > 0 for all z = 1,..., d}. 
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For each j = 1,..., 2*^ we choose for the call, respectively put, profile 
(4.3) rf := —(1 + , respectively rf := 

for some e > 0 for each j = 1,..., 2“^. 

The following proposition enables us to prove analytic dependence of prices of 
basket put and call options on the strike. 

Proposition 4.3. For every j = 1,... ,2"^ let and rf as defined in (14.211 and (|4.3I1 . 
Then for each fc G IR the payoff profile of a basket on a call, respectively put, is 
of the form f^{x) = (e^^ + ... + e^'* — e*^ respectively f^{x) = (e*^ — e^^ — ... — 
e^"^ )^, and 


2d 

(4.4) r\x)=^!l(x) 

i=i 

with fj{x) := f^{x)lQj{x) for a.e. x = {xi,... ,Xd) G R'^. Moreover, for every 
j = 1,...,2'^ and every A: G R we have that x (x) G L^(R'^) and for 

every 2 ; G R'^ the mapping 

(4.5) ke^f^{z — irf) 
has a holomorphic extension. 

Proof. We show the claim for the call option. The put option case is proved anal¬ 
ogously. For every j = 1,...,2'^, under the assumptions of the proposition it is 
elementary to show that x Qpd fj(x) G L^(R'^). Moreover, letting u := z — irj^ 
with 2; G R'^ we have 

JHu)= [ e*<“’*>/^(x)dx 

(4.6) =^k^i{u,k) f gi{n,x} 

JAik) 

with k = {k,... ,k) G R'^ and A{k) := (O^ — fe) n {x G R'^ : /°(x) > 0}. From the 
geometry of A{k) and the holomorphicity of the integrand in ()4.6I) we can deduce 
that the integral in (14.6p is holomorphic in k and thus obtain the assertion of the 
proposition. □ 

4.2 Analyticity Properties of Asset Models 

In this section, we shortly introduce a range of asset models and present analyticity 
properties of their Fourier transforms. For some models and some parameters, the 
domain of analyticity is immediately observable. For some non-trivial cases, we state 
the domain briefly. Throughout the section, T > 0 denotes the time to maturity of 
the option while r > 0 refers to the constant risk-free interest rate. 
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4.2.1 Multivariate Black&Scholes Model 


In the multivariate Black&Scholes model, the characteristic function of LJ, the 
multivariate analogon of (j3.5p . with vr = {b,a) is given by 


(4.7) 




exp 





with covariance matrix a G such that det((T) > 0 and drift b G adhering to 
the drift condition 


(4.8) 



i = 1,... ,d. 


In the Black&Scholes model, analyticity in the parameters is immediately con- 
hrmed, i.e. e- (p^ (z), given in (|4.7I) . is holomorphic for every 2 ; G IR*^. The ad¬ 
missible parameter domain, however, is restricted to parameter constellations such 
that (T is a covariance matrix. 


Remark 4.4. Let rj G be the chosen weight in Conditions \8.1\ and let the open 
set U be given by 

(4.9) U C (0, 00 ) X R X |(T G | a{(j) positive definite^ , 

where a : R‘^('^+l)/2 ^ ]^dxd a{a)ij = 0-(niax{iJ}-l)max{iJ}/2-Hmin{i,j}; 

i,j G for a G R'^('^+^)/2. By construction, cr{(f) is symmetric for any 

a G R'^('^+i)/2. 

Then for every z G R'^, {T,b,5) 1 —>■ pP + given by (14.7p is analytic 

on U. Note that U does not depend on rj. 


4.2.2 Univariate Merton Jump Diffusion Model 

As second Levy model we state the Merton Jump Diffusion model by Merton ( 19761 ). 
The logarithm of the asset price process follows a jump diffusion with volatility 
enriched by jumps arriving at a rate A > 0 with normally M{a, /3^) distributed jump 
sizes. Here, the characteristic function of Llf in (j3.5p with vr = (6, u, a, /3, A) is given 
by 


(4.10) ^gxp 


ibz — 


-^-1 


with drift condition (13.7p turning into 


(4.11) 


a 


b = r---X 




2—1 


where u > 0, a G R and /3 > 0. Since the characteristic function in the Merton 
Jump Diffusion model is composed of analytic functions, it is itself analytic on the 
whole parameter domain. 
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4.2.3 Univariate CGMY Model 


We consider the CGMY model by Carr, Geman, M adan and Yor (2002), which is 


a spec ial case of the extended Koponen’s family of iBovarchenko and Levendorskii 


( 20001 ). with characteristic function (|3.6p of in (13.5p with tt = (6 , C, G, M,Y) 


= exp (T {ibz + Cr(-Y) [(M - izf - + {G + izf - )) , 


where C'>0, G>0, M>0, 0<Y <2 and r(-) denotes the Gamma function. The 
drift 6 E IR is set to 


(4.13) b = r- Cr(-Y) [(M - 1)^ - + [G + Vf - G^] 

to comply with the drift condition (|3.7p for martingale pricing. 

Remark 4.5. For weight r/ E R from Conditions \3fl\ choose an open set U{r]) with 

U(??) C (0,oo) X R X (0,oo) 

X {(G,M) E (0,oo)2 I G-7? > 0, M + r/> 0} X (0,2). 

Then for every z E R, (T, 6, G,G,M,Y) (pP'^={T,b,c,G,M,Y) Jqj. charac¬ 

teristic function ipT of the CGMY model ()4.12p is analytic on U{ri). 

Table 14.21 displays for selected Levy models conditions on the weight ?? E R and 
the index a E (1,2] that guarantee (A3) and the upper bound in (A4) for a fixed 
model parameter constellation. 


Class 

(A3) and the upper bound in 

r] such that 

(A4) hold for 

and a such that 

Brownian Motion 
with drift 


a = 2 

Merton Jump Diffusion 


0 = 2 

Levy jump diffusion with 
characteristics (6, a, F) 

F(dx) < oo 

0 = 2 

univariate GGMY with 
parameters {G,G, M,Y) 
with Y > 1 

r/ E (— min{G, M}, max{G, M}) 

o = Y 


Table 4.2: Conditions on y and a for (A3) and the upper bound in (A4) to hold for a fixed model 
parameter constellation. The selected Levy models are described in more detail in Sections 14.2.IL 


For the GGMY model with parameters C,G,M >0 and Y < 2 it is obvious 
that (|3.11ll is statisfied for ev ery r/ E R with G > rj and M > —rj. Moreover, 
equation (4.12) in Glau ( 20151 ) shows that for Y <1, (|3.12p is satisfied with /3 = 1. 
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4.2.4 Heston Model for Two Assets 


Here we state the two asset version of the multivariate Heston model in the special 
case of having a single, univariate driving volatility process {vt}t>o- The two asset 
price processes are modeled as 


(4.15) 


5/ = and Sf = , for t > 0, 


2„H? 


where H = H^) solves the following system of SDEs, 


dH} = - -cr?J dt + (Tiy/¥tdW^, 

dt + 0-2dlT^^ 
dt>t = k{ 0 - vt) dt + (T^y/m dWf, 

where the Brownian motions Wj, i = 1, 2, 3, are correlated according to W'^) = 
912 , {W^,W^) = /?i 3 , (yV‘^,W^) = P 23 - Following Eberlein, Glau, Papapantoleon (2010), 
the characteristic function of Ht in this framework is 


if 


,P^ = iT,Vo,K,d,ai,o-2,cr3,pi2,pl3,p23) 


(4.16) 


= exp ( Ti 


(z) 


exp 


vq (a — c)(l — exp(—cT)) 


k9 

-2 

^3 L 


(a — c)T — 2 log 


erg 1—5fexp(—cT) 

1 - 5 exp(-cT) 


1-5 


with auxiliary functions 


C = C(^) = - 




/512'7iO'2\ 

ai ) 


z + 


, IZ 


^pi2Cri(T2 

- (afzi + <T2Z2 + 2pi2fTifT2Zi2;2 + icrfzi + io-2Z2) , 

(4.17) « = = K- ipuO-iasZi - ip230-20-3Z2, 

c = c(z) = ^Ja{zf -alciz), 

a{z) - c{z) 

9 = 9{z) = . ^ . X , 

a{z) + c{z) 

and positive parameters vq, k, 9 and fulfilling the Feller condition 


(4.18) 


(t| < 2k9, 


ensuring an almost surely non-negative volatility process {vt)t>o- Obviously, for 
each z G IR,'^, the characteristic function ipP^=('^’‘^o,K,e,ai,a2,a3,pi2,pi3,p23)(^^'j Qf (| 4 ,i 6 p 
is analytic in vq and 9. Let us additional l y men tion that the analiticity in z has 
already been investigated by Levendorskii ( 20121 . Section 2.3, Lemma 2.1) for the 
univariate Heston model. 
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5 Numerical Experiments 


We apply the Chebyshev interpolation method to parametric option pricing con¬ 
sidering a variety of option types in different well known option pricing models. 
Moreover, we conduct both an error analysis as well as a convergence study. The 
first focuses on the accuracy that can be achieved with a reasonable number of 
Chebyshev interpolation points. The latter confirms the theoretical order of conver¬ 
gence derived in Section O when the number of Chebyshev points increases. Finally, 
we study the gain in efficiency for selected multivariate options. 

We measure the numerical accuracy of the Chebyshev method by comparing 
derived prices with prices coming from a reference method. We employ the reference 
method not only for computing reference prices but also for computing prices at 
Chebyshev nodes Price^ ’ ’ with {ki,... jko) £ J during the precomputation 
phase of the Chebyshev coefficients Cj, j G J, in (|2.9I) . Thereby, a comparability 
between Chebyshev prices and reference prices is maintained. 

We implemented the Chebyshev method for applications with two parameters. 
To that extent we pick two free parameters pi^, pi^ out of (13.2p . I < h < 12 < D, 
in each model setup and fix all other parameters at reasonable constant values. 
We then evaluate option prices for different products on a discrete parameter grid 
^ ^ X defined by 

^ = { ^) , hi, h2 G {0,... , 100}I , 

(5-1) k- k- / ^ 

Pi-' = Pi, + [Pii -EiJ > e • • • > 100}> j e 2 }. 

Once the prices have been derived on V, we compute the discrete L°°{V) and L‘^{V) 
error measures. 


eL°°{N) = m^ Price^ — Ij;^{Priceh){p) 
p&V 


(5.2) 


ei2(iV) = 


\ 


Ap ^ jphce^ — I-^{Price^'^){p) 


p&V 


where = *100 ^ "^100 ^ ; fo interpret the accuracy of our implementation and 

of the Chebyshev method as such. 


5.1 European Options 

We consider a plain vanilla European call option on one asset as well as a European 
digital down&out option, first. Both derivatives have been introduced in Table HTl 
For these products we investigate the performance of the Chebyshev interpolation 
method for the Heston model and the Levy models of Black&Scholes, Merton and the 
CGMY model. We keep the strike parameter constant, p^ = k = \og{K) for K = 1. 
As previously discussed in Remark 13.81 this is no restriction of generality. We also 
disregard interest rates, setting r = 0. For the three Levy models we vary the matu¬ 
rity T (in years) as well as the option moneyness Sq/K. All three models fall within 
the scope of Corollary 13.81 where the error is analyzed explicitly. Thus we expect 
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Model 

fixed parameters 


free 

parameters 


P 



pZ 


pi 



p2 

BS 

K = 

1 

a = 

0.2 

So/K 

G [0.8, 

1.2] 

T 

G [0.5,2] 

Merton 

K = 

1 

a = 

0.15, 

So/K 

G [0.8, 

1.2] 

T 

G [0.5,2] 




a = 

-0.04, 









/3 = 

0.02, 









A = 

3 






CGMY 

K = 

1 

C = 

: 0.6, 

So/K 

G [0.8, 

1.2] 

T 

G [0.5,2] 




G = 

:10, 









M -- 

= 28, 









Y = 

: 1.1 






Heston 

K = 

1 

T = 

2, 

So/K 

G [0.8, 

1.2] 

Vo 

G [0.l2, 0 . 42 ] 




K = 

1.5, 









e = 

0.22, 









a = 

0.25, 









P = 

0.1 







Table 5.1: Parametrization of models and European call option. 


(sub)exponential convergence for the three Levy models. For the Heston model we 
vary Sq/K and let vq as one of the model parameters float. Due to the analyticity 
of the Fourier transform of the payoff in Sq/K and of the characteristic function of 
the process in vq, compare Section [4.2.41 we expect this convergence for the Heston 
model as well. A detailed overview of the realistically chosen parametrization is 
given by Table 15.11 For numerical integration in Fourier pricing we use Matlab’s 
quadgk routine over the interval [0,oo) with absolute precision bound of e < 

The first question we address concerns the achievable accuracy with a fixed 
number of Chebyshev polynomials. We set Ni = N 2 = 10 and precompute the 
Chebyshev coefficients as defined in ()2.9p with D = 2 using the parametrization 
of Table 15.11 for the models therein. We evaluate the resulting polynomial over a 
parameter grid of dimension D = 2 and compute the approximate European option 
prices in each node. As a comparison, we also compute the respective Fourier price 
via numerical integration of the accordingly parametrized integrand in (I3.4p . Fig¬ 
ure 0 shows the results for the European call option. The accuracy achieved by 
N = Ni = N 2 = 10 shows a significant spread over the four different models but 
reaches very satisfying error levels from 10“^ to 10“^*^. Increasing the number of 
Chebyshev points further improves the accuracy. Since at its core the implemen¬ 
tation of the Chebyshev method consists of summing up matrices, this refinement 
comes at virtually no additional cost. 

In the same parametrization setting we price a European digital down&out op¬ 
tion. While a call payoff profile is not differentiable but at least continuous, the 
digital payoff function is not even continuous, compare Table 14.11 This reduc¬ 
tion in smoothness of the payoff function reduces the accuracy of the interpola¬ 
tion p I—>■ Price^P^ as well. We compare the Chebyshev interpolation again with 
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„ Chebyshev price error, BS 

xlO 


5 ^ 



„ Chebyshev price error, CGMY 

xlO 


5 





Figure 5.1: Absolute pricing error for a European call option with strike K = 1 in various models. 
We compare the Chebyshev interpolation with N = Ni = = iQ to classic Fourier pricing by 

numerical integration. The parametrization of the models and the option has been chosen according 
to Table [5T] We observe that the achieved accuracy varies significantly between the models. The 
order of accuracy of Black&Scholes prices is identical to the accuracy for the CGMY case. The 
Merton model falls behind by one order of magnitude, while the Heston model shows a very strong 
fit between prices and their approximation by the Chebyshev method. 


Ni = N 2 = 10 to classic Fourier pricing by numerical integration. The parametriza¬ 
tion of the models and the option has been chosen again according to Table 15.11 
where K = 1 now denotes the parameter of the digital option. Figure shows the 
results. Comparing the results of the call option pricing in Figure EH to the results 
of the digital option pricing in Figure 15.21 we observe that the accuracy achieved 
by A^i = N 2 = 10 is reduced by a factor 10^ to 10^. This demonstrates how little 
the reduced smoothness of the payoff profile effects the accuracy of the Chebyshev 
method in these cases. Furthermore we conduct an empirical convergence study for 
this very same setting of option and model parametrization. For an increasing de¬ 
gree N = Ni = N 2 , the Chebyshev polynomial is set up and prices over a parameter 
grid of structure (j5.19ll are computed. Again, Fourier pricing serves as a comparison. 
For each N G {1,...,30}, the error measures £ l °° and £^ 2 , defined by (|5.2p on the 
discrete parameter grid V defined in (15.1911 . are evaluated. We observe exponential 
convergence for all four models. Figure 15.31 shows the decay for the European call 
option while Figure EH does so for the European digital down&out option. 


Defining and setting = C + \/C^ ~ 1; the theoretical convergence analysis 

predicts a slope of the error decays in both Eigure 15.31 as well as Eigure 15.41 of at 
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Chebyshev price error, Merton 
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Figure 5.2: Absolute pricing error for a European digital down&out option with A = 1 in various 
models. The plots of this figure correspond to the plots of Figure [5TT] where a Europan call option 
was priced instead. We compare again the performance of the Chebyshev method for = A^i = 
N 2 = 10. All four price surfaces lose between one to two orders of magnitude in accuracy compared 
to their call option counterparts. 



Chebyshev N 

Figure 5.3: Convergence study for the Black&Scholes, Merton, CGMY and the Heston model for 
prices of a European call option parametrized as stated in Table [57T1 The reference price is derived 
by Eourier pricing and numerical integration with an absolute accuracy of 10“^“^, which is reached 
by all models for N = Ni = N 2 « 25 the latest. The error decay for the three Levy models of 
Black&Scholes, Merton and the CGMY model roughly coincide, extending the hndings from Figure 
[5T] 
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Chebyshev error decay 



Chebyshev N 

Figure 5.4: Convergence study for prices of a European digital down&out option in the 
Black&Scholes, Merton, CGMY and the Heston model parametrized as stated in Table 15.11 The 
results of this figure correspond to the results shown in Figure [5^ where the error decay of prices 
of a European call option in the same model setup was analyzed. Now higher N and thus more 
Chebyshev nodes are needed to reach the same levels of accuracy as before. 


least 

•s = logio ie) ~ -0.48 

or steeper. Note that G = leads to a higher p and for this analysis we 

boK—boK 

compare the slope to the minimal g value as in Corollarv l2.3l Empirically, we observe 
a slope for the Black&Scholes model of about 5 bs = —0.64, for the Merton model 
of ^Merton = —0.61 and for the CGMY model of 5cgmy = —0.61. Thus, the error 
in each Levy model empirically confirms the theoretical claim of Remark 13.81 


5.2 Basket and Path-dependent Options 


In this section we use the Chebyshev method to price basket and path-dependent 
options. First, we apply the method to interpolate Monte-Carlo estimates of prices 
of financial products and check the resulting accuracy. To this aim we exemplar- 
ily choose basket, barrier and lookback options in 5-dimensional Black&Scholes, 
Heston and Merton models. Second, we combine the Chebyshev method with a 
Crank-Nicolson finite differen ce solver with Brennan Schwartz approximation, see 
Brennan and Schwartz! ()l977l L for pricing a univariate American put option in the 
Black&Scholes model. 

In our Monte-Carlo simulation we use 10® sample paths, antithetic variates as 
variance reduction technique and 400 time steps per year. The error of the Monte- 
Carlo method cannot be computed directly. We thus turn to statistical error analysis 
and use the well-known 95% confidence bounds to determine the accuracy. Following 
the assumption of a normally distributed Monte-Carlo estimator with mean equal 
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to the estimator’s value and variance equal to the empirical variance of the payoff 
on the Monte-Carlo samples these bounds are derived. The confidence bounds then 
yield a range around the mean that includes the true price with 95% probability. 
We pick two free parameters pq, pi^ out of (13.211 . 1 < ii < i 2 < D, 'm each model 
setup and fix all other parameters at reasonable constant values. In this section we 
define the discrete parameter grid V C x \p^^-,Pi^ by 

{ (Pil'^pt^) > 40} I , 

^ ^ {0,..., 40}, j E {1, 2}, 

and call V test grid. On this test grid the largest confidence bound is 0.025 an 
on average lees than 0.013. For the finite difference method we find that the abso¬ 
lute error between numerical approximation and option price is below 0.005 on all 
computed parameter tuples in V. This error bound was computed by comparing 
each approximation to the limit of the sequence of finite difference approximations 
with increasing grid size. In our calculations we work with a grid size in time as 
well as in space (log-moneyness) of 50 • max{l,T} and compared the result to the 
resulting prices using grid sizes of 1000 • max{l,T}. This grid size has been deter¬ 
mined as sufficient for the limit due to hardly any changes compared to grid sizes of 
500 • max{l, T}. 

Here, our main concern is the accuracy of the Chebyshev interpolation when we 
vary for each option the parameters strike and maturity analogously to the previous 
section. For N E {5,10,30}, we precompute the Chebyshev coefficients as defined 
in (12.9p with D = 2 where always Ni = N 2 = N. An overview of fixed and free 
parameters in our model selection is given in Table[52J For computational simplicity 
in the Monte-Carlo simulation, we assume uncorrelated underlyings. 

Let us briefly define the payoffs of the multivariate basket and path-dependent 
options. The payoff profile of a basket option for d underlyings is given as 


(5.3) 


V = 

ki. 

Pi/ = 





We denote St = {St,...,Sf), 5^ := mino<t<T <5/ and S/ := msiXQ<t<TS{. A 
lookback option for d underlyings is defined as 


/'f(4....,4) = -A'j . 

As a multivariate barrier option on d underlyings we define the payoff 

/'^ ({S(t)}„<,<T) = f f 1E j - if j ■ ■ 

For an American put option the payoff is the same as for a European put. 


f^{St)={K-St)+, 


30 


Model 

fixed parameters 

1 2 

P P 

free parameters 

1 2 

P P 

BS 

sl = 100, 

r = 0.005 


= 0.2 

K G [83.33, 125] 

T G [0.5,2] 

Heston 

= 100, 
r = 0.005 

Kj — 2, 

e.j = 0.22, 
aj = 0.3, 

Pj = -0.5, 

Vjfi = 0.22 

K G [83.33, 125] 

T G [0.5,2] 

Merton 

si = 100, 

r = 0.005 

^3 

aj 

^3 

^3 

= 0.2, 

= -0.1, 

= 0.45, 

= 0.1 

K G [83.33, 125] 

T G [0.5,2] 


Table 5.2: Parametrization of models, basket and path-dependent options. The model parameters 
are given for j = 1,..., d to reflect the multivariate setting with free parameters strike K and 
maturity T. Note that in contrast to the two dimensional Heston model described in Section f4. 2. 41 
we use here in the numerical experiments a multivariate Heston model in which the volatility of 
each underlying is driven by its own volatility process. 


but the option holder has the right to exercise the option at any time t up to matu¬ 
rity T. 


Model 

Option 

£1,00 

MC price 

MC conf 

bound 

Cl price 

BS 

Basket 

1.338 • 

10-^ 

8.6073 

1.171 • 

10-2 

8.4735 

Heston 

Basket 

9.238 • 

10-2 

0.0009 

1.036 • 

10-^ 

0.0933 

Merton 

Basket 

9.815 • 

10-2 

8.8491 

1.552 . 

10-2 

8.7510 

BS 

Lookback 

2.409 • 

10-^ 

9.4623 

9.861 • 

10-^ 

9.2213 

Heston 

Lookback 

5.134 . 

10-1 

0.0314 

6.472 . 

10-^ 

-0.4820 

Merton 

Lookback 

2.074 . 

10-1 

1.0919 

9.568 • 

10-3 

0.8844 

BS 

Barrier 

1.299 • 

10-1 

1.0587 

5.092 . 

10-3 

1.1887 

Heston 

Barrier 

1.073 • 

10-1 

2.7670 

9.137- 

10-3 

2.6597 

Merton 

Barrier 

9.916 • 

10-2 

1.3810 

1.102 . 

10-2 

1.4802 


Table 5.3: Interpolation of exotic options with Chebyshev interpolation. N = 5 and d = 5 in all 
cases. In addition to the L°° errors the table displays the Monte-Carlo (MC) prices, the Monte- 
Carlo confidence bounds and the Chebyshev Interpolation (Cl) prices for those parameters at which 
the L°° error is realized. 


We now turn to the results of our numerical experiments. In order to evaluate 
the accuracy of the Chebyshev interpolation we look for the worst case error sloo. 
The absolute error of the Chebychev interpolation method can be directly computed 
by comparing the interpolated option prices with those obtained by the reference 
numerical algorithm i.e. either the Monte-Carlo or the Finite Difference method. 
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Model 

Option 


MC price 

MC conf 

bound 

Cl price 

BS 

Basket 

2.368 • 

10"^ 

2.4543 

7.493 • 

10-3 

2.4566 

Heston 

Basket 

2.134- 

10"3 

3.1946 

1.073 • 

10-2 

3.1925 

Merton 

Basket 

3.521 • 

10"3 

6.1929 

2.231 • 

10-2 

6.1894 

BS 

Lookback 

2.861 • 

10"^ 

0.9827 

4.197- 

10-3 

0.9541 

Heston 

Lookback 

1.098 • 

10-1 

2.0559 

4.826 - 

10-3 

2.1656 

Merton 

Lookback 

3.221 • 

10-2 

4.7072 

1.264 - 

10-2 

4.7394 

BS 

Barrier 

4.414 • 

10-^ 

5.3173 

1.725 - 

10-2 

5.3129 

Heston 

Barrier 

5.393 • 

10-3 

0.7158 

5.879 - 

10-3 

0.7212 

Merton 

Barrier 

3.376 • 

10-3 

9.2688 

2.302 - 

10-2 

9.2722 


Table 5.4: Interpolation of exotic options with Chebyshev interpolation. = 10 and d = 5 in all 
cases. In addition to the L°“ errors the table displays the Monte-Carlo (MC) prices, the Monte- 
Carlo confidence bounds and the Chebyshev Interpolation (Cl) prices for those parameters at which 
the L°° error is realized. 


Model 

Option 

£L°° 

MC price 

MC conf 

bound 

Cl price 

BS 

Basket 

1.452 - 

10-3 

5.1149 

1.200 - 

10-2 

5.1163 

Heston 

Basket 

1.047- 

10-3 

7.6555 

1.371 - 

10-2 

7.6545 

Merton 

Basket 

3.765 - 

10-3 

7.2449 

2.359 - 

10-2 

7.2412 

BS 

Lookback 

3.766 - 

10-3 

25.9007 

1.032 - 

10-2 

25.9045 

Heston 

Lookback 

1.914 - 

10-3 

16.4972 

9.754 - 

10-3 

16.4991 

Merton 

Lookback 

3.646 - 

10-3 

27.1018 

1.623 - 

10-2 

27.1054 

BS 

Barrier 

5.331 - 

10-3 

5.6029 

1.730 - 

10-2 

5.6082 

Heston 

Barrier 

2.486 - 

10-3 

3.6997 

1.353 - 

10-2 

3.6972 

Merton 

Barrier 

4.298 - 

10-3 

6.6358 

2.309 - 

10-2 

6.6315 


Table 5.5: Interpolation of exotic options with Chebyshev interpolation. = 30 and d = 5 in all 
cases. In addition to the L°° errors the table displays the Monte-Carlo (MC) prices, the Monte- 
Carlo confidence bounds and the Chebyshev Interpolation (Cl) prices for those parameters at which 
the L°° error is realized. 


Since the Chebychev interpolation matches the reference method on the Chebychev 
nodes, we will use the out-of-sample test grid as in (|5.3D . Table 15.31 shows the 
numerical results for the basket and path-dependent options for N = 5, Table 15.41 
for = 10 and Table 15.51 for N = 30. In addition to the errors the tables 
display the Monte-Carlo (MC) prices, the Monte-Carlo confidence bounds and the 
Chebyshev Interpolation (Cl) prices for those parameters at which the L°° error is 
realized. 

The results show that for = 30 the accuracy is for all selected options at 
a level of 10“^. We see that the Chebyshev interpolation error is dominated by 
the Monte-Carlo confidence bounds to a degree which renders it negligible in a 
comparison between the two. For basket and barrier options the L°° error already 
reaches satisfying levels of order 10“^ at A^ = 10 already. Again, the Chebyshev 
approximation falls within the confidence bounds of the Monte-Carlo approximation. 
Thus, Chebyshev interpolation with only 121 = (10-|-1)^ nodes suffices for mimicking 
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the Monte Carlo pricing results. This statement does not hold for lookback options, 
where the L°° error still differs noticeably when comparing = 10 to iV = 30. 
As can be seen from Table 15.31 Chebyshev interpolation with = 5 may yield 
unreliable pricing results. For lookback options in the Heston model we even observe 
negative prices in individual cases. Chebyshev pricing of American options in the 


N 

eL°° 

FD price 

Cl price 

5 

3.731 • 10-^ 

1.9261 

1.9224 

10 

1.636 • 10-3 

12.0730 

12.0746 

30 

3.075 • 10-3 

6.3317 

6.3286 


Table 5.6: Interpolation of one-dimensional American puts with Chebyshev interpolation in the 
Black&Scholes model. In addition to the errors the table displays the Finite Differences (FD) 
prices and the Chebyshev Interpolation (Cl) prices for those parameters at which the L°° error is 
realized. 


Black&Scholes model is even more accurate as illustrated in Table (5^ Here, already 
for N = 5 the accuracy of the reference method is achieved. We conclude that the 
Chebyshev interpolation is highly promising for the evaluation of multivariate basket 
and path-dependent options. Yet the accuracy of the interpolation critically depends 
on the accuracy of the reference method at the nodal points which motivates further 
analysis that we perform in the subsequent subsection. 


5.2.1 Interaction of Approximation Errors at Nodal Points and Inter¬ 
polation Errors 

The Chebyshev method is most promising for use cases, where computationally in¬ 
tensive pricing methods are required. Then, for computing the prices at the Cheby¬ 
shev nodes in order to set up the interpolation, the issue of distorted prices at the 
Chebyshev nodes and their consequences rises naturally. The observed noisy prices 
at the Chebyshev nodes are 

Pricef^ . 

where jg approximation error introduced by the underlying numerical 

technique at the Chebyshev nodes. Due to linearity, the resulting interpolation is of 
the form 


(5.4) 

/^(PnceP)(p) 

= lN{Pricek^){v) + 

with the 

error function 


(5.5) 

No 

jD=0 

Ni 

■ ■ ■ 

ji=0 

with the coefficients for j = (ji 

,...,jD) ^ J given by 

(5.6) 


Ni No D 

L L li' 

fci=0 k]j=0 2=1 

^ Ui Ni ) 

2=1 


ki 
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If < £ fQp all Chebyshev nodes -^g obtain 


(5.7) 


D 

e(p)| <2^eJ](iVi + l), 

i=l 


since the Chebyshev polynomials are bounded by 1. This yields the following remark. 




Remark 5.1. Let V 3 p Price^ be given as in Theorem \2.S\ and assume that 
< £ for all Chebyshev nodes _ Then 


m^ |Pnce^ — Ij;^{Priee^^^){p)\ 


(5.8) 


D 


D 


D 


< 2^+^ • V ■ 


E a n 7—^ + +!)■ 

i=i j=i / i=i 


The following example shall illustrate the practical consequences of Remark 15.11 
In the setting of Corollary 1,1.81 we set \So / K, So/ K] = [0.8,1.2], [T, T] = [0.5,2]. 
This results in = | and = -^ = 5. Thus, for gi = 2.9 G (1,3) and 

g 2 = 9.8 G (1, 5 + \/^), Remark l5 .1 1 yields with A^i = A ^2 = 6, 

max |Przce^ — Ij^{Price^''l){p)\ < 0.0072 + 196 • e. 


In this example, the accuracy of the reference method has to reach a level of 10“® 
to guarantee an overall error of order 10“^. This demonstrates a trade-off between 
increasing A^i and N 2 compared to the accuracy of the reference method. The error 
bound above is rather conservative. Our experiments from the previous section sug¬ 
gest that this bound highly overestimates the errors empirically observed. However, 
the presented error bound from Remark 15.11 can guarantee a desired accuracy by de¬ 
termining an adequate number of Chebyshev nodes and the corresponding accuracy 
of the reference method used at the Chebyshev nodes. For practical implementation 
we suggest the following procedure. For a prescribed accuracy the Ni, i = 1,..., D, 
can be determined from the first term in (15.8p by choosing Ni, i = 1,..., D, as small 
as possible such that the prescribed accuracy is attained. Accordingly, the accuracy 
that the reference method needs to achieve is bounded by the second term. A very 
accurate reference method in combination with small Ni, i = 1,... ,D, promise best 
results. With this rule of thumb in mind the experiments of Section 15.3.21 below 
have been conducted. 


5.3 Study of the gain of efficiency 

In the previous section we investigated the accuracy of the Chebyshev polynomial 
interpolation method using Fourier, Monte-Carlo and Finite Difference as reference 
pricing methods. Finally, we investigate the gain in efficiency achieved by the method 
in comparison with Fourier pricing as well as in comparison to Monte Carlo pricing. 
In Section 15.3.11 we compute the results on a standard PC with an Intel i5 CPU, 
2.50 GHz with cache size 3 MB. In Section 15. 3.2 1 we used a PC with Intel Xeon CPU 
with 3.10 GHz with 20 MB SmartCache. All codes are written in Matlab R2014a. 
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5.3.1 Comparison to Fourier pricing 

Here, we compare the method to Fourier pricing. We choose the pricing problem of 
a call option on the minimum of two assets as an example. Today’s values of the 
underlying two assets are fixed at 

(5.9) = 1, = 1.2. 

Modeling the future development of the underlyings, {Sl)t>o, j G {1,2}, we con¬ 
sider two bivariate models, separately. First, the two assets will be driven by the 
bivariate Black&Scholes model of Section S2TJ The bivariate Black&Scholes model 
is parametrized by a covariance matrix a G that we choose to be given by 

(5.10) an = 0.2^, a^ = 0.01, 0-22 = 0.25^. 

In a second efficiency study, asset movements follow the more involved bivariate Hes¬ 
ton model in the version of Section [4. 2. 41 above for which we choose the parametriza- 
tion 


(5.11) 


Vo = 0.05, 

K = 0.4963, 
0 = 0.2286, 


ai = 0.15, 

(T2 = 0.2, 
0-3 = 0.1, 


Pi3 = 0.01, 
P 12 = 0, 

P 23 = 0.02. 


In both cases we neglect interest rates, thus setting r = 0. The benchmark method, 
that is Fourier pricing, is evaluated using Matlab’s quad2d routine. We prescribe 
an absolute and relative accuracy of at least 10“® from the integration result and 
integrate the Fourier integrand over the domain = [—50, 50] x [0,50], with a maxi¬ 
mum number of 4000 function evaluations. We use Fourier integration with the same 
accuracy specihcations at the Chebyshev nodes to set up the Chebyshev method for 
pricing based on strike K and maturity T as the two free parameters taking values 
in the intervals 


(5.12) 


^ G [.fiimin, 
T G [Tjninj 



Krain = 0.8, ifimax = 1-2, 
Train = 0.5, Tmax = 2. 


For a fair comparison, the number of Chebyshev polynomials is chosen such that 
Chebyshev interpolation prices yield an accuracy that matches the accuracy of the 
benchmark method resulting in 

(5.13) lV^h®eby = ll and = 23, 

for the bivariate Black&Scholes model and the bivariate Heston model, respectively. 
Figure [531 illustrates the absolute errors over the whole K xT domain of inter¬ 
est between Fourier pricing and the Chebyshev method for both models, with the 
Chebyshev interpolator being based on -|-1 polynomials in the Black&Scholes 

model case and + 1 polynomials in the Heston model case. 

In order to set up an efficiency study, we will compute sets of prices with in¬ 
creasing number of parameter tupels. To this end, when the offline phase of the 
Chebyshev method has been completed we compute 98 pricing surfaces, that is for 
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xlO 

2 


Chebyshev accuracy, BS 


xlO’“ 

5 


Chebyshev accuracy, Heston 



Figure 5.5: Left: Difference between prices from the Fourier method and Chebyshev interpolation 
in the bivariate Black&Scholes model over the whole parameter domain of interest. The model 
is parametrized as indicated by (15.1011 . Chebyshev interpolation is based on A^cheby + 1 = 12 
Chebyshev polynomials. Right: The respective plot for the Heston model parametrized as in (IstTIi . 
Here, Chebyshev interpolation is based on A^cSby" + 1 = 24 Chebyshev polynomials. We achieve 
an absolute accuracy of order 10“® in both cases, thus matching the accuracy that the benchmark 
method Fourier pricing reaches. 


each M G {3,..., 100} we compute prices for all parameter tuples from Qm defined 
by 

0M = {(i^f ) I 

(5.14) Tf = Tmin + ^ (Tmax “ T^in), foT 1 <i,j < M}. 

The computation time consumed by the Chebyshev offline phase is stored. Also, 
for each M G {3,..., 100}, run-times for deriving all |0 m| = prices are measured 
and stored for both routines, the Fourier pricing method and the Chebyshev inter¬ 
polation algorithm. Figure 15.61 depicts these run-time measurements and Table 15.71 
provides a second perspective. 

In the Black&Scholes model case, the offline phase required = 8 seconds 

for deriving option prices at all (A^cheby d" = 14,4 Chebyshev nodes. The Heston 
model required = 101 seconds for the (A^cheby" + = 576 supporting 

prices. Taking this initial investment into account deems pricing with the Chebyshev 
method rather costly when only few option prices are derived after the offline phase 
has been completed. Yet, as Figure (5^ shows, the increase in pricing speed that is 
achieved once the Chebyshev algorithm has been set up eventually outpaces Fourier 
pricing. From our experiments we conclude that the Chebyshev method already 
outruns its benchmark Fourier pricing in terms of total run-times when the number 
of prices to be computed exceeds (A'cheby + or (-^cSby” + respectively. 
Additionally, Table 15.71 highlights that in both models for a total number of 50^ 
parameter tuples, the Chebyshev method exhibits a significant decrease in (total) 
pricing run-times. For the maximal number investigated, i.e. 100^ parameter tuples, 
pricing in the Black&Scholes model results in 95% of run-time savings and 90% run¬ 
time savings in the Heston model case in our implementation. This results from the 
fact the online phase in the Chebyshev method consist of computationally cheap 
polynomial evaluations and elementary assembling. 
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Black&Scholes Call option: Offline + online time 



5 15 25 35 45 55 65 75 85 95 

M (evalnation of NP different parameter pairs) 


Heston Call option: Offline + online time 



5 15 25 35 45 55 65 75 85 95 

A/ (evaluation of KP different parameter pairs) 


Figure 5.6: Comparison of pricing times between Fourier pricing and the Chebyshev method for a 
call option on the minimum of two assets in the Black&Scholes model (left) and the Heston model 
(right). For each M € {3,..., 100}, run-times for deriving option prices for all parameter 
tupels from 0 m defined by (I5d4ll are depicted. In both model cases, computation times for the 
Chebyshev method contain the duration of the offline phase that has to be conducted once in the 
beginning. The Fourier and the Chebyshev curves roughly intersect when M = -|- 1 = 12 for 

the Black&Scholes model and when M = -|- 1 = 24 for the Heston model, respectively. 


BS Heston 


M 

10 

50 

75 

100 

10 

50 

75 

100 

rpChehy / n 
-^ online 

0.18 

4.54 

10.20 

18.11 

0.70 

17.58 

39.66 

69.82 

rpChehy / \ 

offline+online 

8.06 

12.42 

18.07 

25.98 

101.96 

118.85 

140.92 

171.08 

2^Fourier 

5.34 

131.96 

301.82 

528.74 

17.60 

442.62 

991.33 

1788.08 

ry-iCheby 
^ offline+online 
2^Fourier 

151% 

9.41% 

5.99% 

4.91% 

579.27% 

26.85% 

14.22% 

9.57% 


Table 5.7: Efficiency study for the bivariate Black&Scholes model and the bivariate Heston model: 
A selection of the results fully depicted in Figure 1531 With increasing number of computed prices, 
the Chebyshev algorithm increasingly profits from the initial investment of the offline phase. 
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5.3.2 Comparison to Monte-Carlo pricing 

In this section we choose a multivariate lookback option in the Heston model, based 
on 5 underlyings, as example. For the efficiency study we first vary one parameter, 
then we vary two. 

Variation of one model parameter 

For the multivariate lookback option in the Heston model the following parameters 
are fixed with j = 1,..., 5 as 

si = 100, r = 0.005, K = 100, T = 1, 

(5 15) u ’ ; ; ’ 

kj = 2, 0j=O. 2^, pj =-0.5, vj,o = 0.2^. 

As free parameter in the Chebyshev interpolation we pick the volatility of the volatil¬ 
ity coefficient a = (Tj, j = 1,..., 5, 


(5.16) <7 € [cTminj <7max]) <^min — 0.1, i7jnax — 0.5. 

The benchmark method is the Monte-Carlo pricing, again with 10® sample paths, 
antithetic variates as variance reduction technique and 400 time steps per year. We 
refer to this setting as benchmark setting. 

Following the discussion from Section 15.2.H for the evaluation of the prices at 
the nodal points we guarantee a small e by the Monte-Carlo method we enrich 
the Monte-Carlo setting to 5 • 10® sample paths, antithetic variates and 400 time 
steps per year. In Table 15.81 we present the accuracy results for the Chebyshev 
interpolation with A^cheby^ ~ ® based on the enriched Monte-Carlo setting. To this 
aim, we compare the absolute differences of the Chebyshev prices and the enriched 
Monte-Carlo prices on the test grid P C [p,p]. 


(5.17) 


P= A:E{0,...,20}}, 

k 

^ — <7min T 7^ (<7max <7min) j k E {0, . . . , 20}. 


The highest observed error on the test grid is at a level of 10“^. On the same 
test grid the benchmark Monte-Carlo setting has a worst case confidence bound of 
1.644 • 10“^ and by comparing the benchmark Monte-Carlo prices to the enriched 
Monte-Carlo prices on this test grid, the maximal absolute error is 7.361 • 10“^. 
Therefore, we conclude that the Monte-Carlo benchmark setting and the presented 
Chebyshev interpolation method have a roughly comparable accuracy and on the 
basis of this accuracy study we now turn to the comparison of run-times. 

We compare run-times of the Chebyshev interpolation with = 6, in which 

the offline phase is based on the enriched Monte-Carlo setting, to the run-times of 
the Monte-Carlo benchmark setting described above. 

Table 15.91 provides the respective results. The results for M = 1 have been 
empirically measured, all others have been extrapolated from that since for each 
parameter set, the same amount of computation time would have to be invested. 
The table indicates that from M = 50 onwards interpolation by Chebyshev is faster. 
In Figure 15.71 we present additionally for each M = 1,..., 100 the run-times of 
the Chebyshev interpolation method, including the offline phase, compared to the 
Monte-Carlo method. Here we observe that for M = 35 both lines intersect and for 
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Varying 

£L^ 

MC price 

MC conf. bound 

Cl price 

a 

9.970-10“® 

18.6607 

4.592 • 10“® 

18.6707 


Table 5.8: Interpolation of multivariate lookback options with Chebyshev interpolation for N = 6 
based on an enriched Monte-Carlo setting with 5 • 10® sample paths, antithetic variates and 400 
time steps per year. In addition to the L°° error on the test grid we also report the Monte-Carlo 
(MC) price, the Monte-Carlo confidence bound and the Chebyshev Interpolation (Cl) price for 
those parameters at which the L°° error is realized. We observe that the accuracy of the Chebyshev 
interpolation for = 6 is roughly in the same range as the accuracy of the benchmark Monte-Carlo 
setting (worst case confidence bound of 1.644 • 10~^ and worst case error of 7.361 • 10“®). 


M > 35 the Chebyshev interpolation method is faster. Contrary to Section 15.3.11 
the intersection of the two lines does not occur at M = -|-1. This reflects the 

fact that in the offline phase for the Chebyshev interpolation a Monte-Carlo method 
with more sample paths has been used. 


Heston Lookback option: Offline -f- online time 



Figure 5.7: EfSency study for a multidimensional lookback option in the Heston model with 5 
underlyings varying one model parameter a. Comparison of run-times of Monte-Carlo pricing with 
Chebyshev pricing including the offline phase. Both methods have been set up to deliver comparable 
accuracies. We observe that both curves roughly intersect at M = 35. 


Variation of two model parameters 

Now, we vary two parameters, we choose pj = p, j = 1,..., 5, and vary 


(5.18) 


P £ [Pmin) Pmax] j Pmin — 1 j Pmax — 1 j 

T G [(Tmin, CTmax]) ^min — 0.1, Tmax — 0.5, 


fixing all other parameters to the values of setting (j5.15l) . Guaranteeing a roughly 
comparable accuracy between the Chebyshev interpolation method and the bench¬ 
mark Monte-Carlo pricing, we use the following test grid V C [fTmim Tmax] x [Pmim /O ma x], 

r= A:i,A:2G{0,...,20}}, 

A;i 

(5.19) (T ® = (Tniin T (<^max Tmin) , k\ G {0, . . . , 20}, 

^2 

P ^ ~ Pmin T (Pmax Pmin) ; ^2 G {0, . . . , 20}. 
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M 

1 

10 

50 

100 

^Cheby / x 
^online 

2.7-10“5 

2.7-10“^ 

1.4-10“3 

2.7-10“3 

T^Cheby / x 

^ offline+online ' 

1.2-10^ 

1.2-104 

1.2-104 

1.2-104 

^Monte-Carlo 

3.4-102 

3.4-10^ 

1.7-104 

3.4-104 

^-iCheby 

offline-honline 

^Monte-Carlo 

3473.4% 

347.3% 

69.5% 

34.73% 


Table 5.9: Efficiency study for a multivariate lookback option in the Heston model based on 5 
underlyings. Here, we vary one model parameter and compare the Chebyshev results to Monte- 
Carlo. Both methods have been set up to deliver comparable accuracies. With increasing number 
of computed prices, the Chebyshev algorithm increasingly profits from the initial investment of the 
offline phase. 


Varying 


MC price 

MC conf. bound 

Cl price 

(T, p 

5.260 - 10-^ 

5.239 

1.428 - 10-^ 

5.292 


Table 5.10: Interpolation of multivariate lookback options with Chebyshev interpolation for N = 6 
based on an enriched Monte-Carlo setting with 5 ■ 10® sample paths, antithetic variates and 400 
time steps per year. In addition to the Z/°° error on the test grid we also report the Monte-Carlo 
(MC) price, the Monte-Carlo confidence bound and the Chebyshev Interpolation (Cl) price for 
those parameters at which the L°° error is realized. We observe that the accuracy of the Chebyshev 
interpolation = 6 is roughly in the same range as the accuracy of the benchmark Monte-Carlo 
setting (worst case confidence bound of 6.783 • 10“^ and worst case error of 2.791 • 10“^). 


In Table 15.101 we present the accuracy results for the Chebyshev interpolation with 
-^Cheby'^ = 6 based on the enriched Monte-Carlo setting. Comparing the benchmark 
Monte-Carlo setting and the enriched Monte-Carlo setting on this test grid, we 
observe that the maximal absolute error is 2.791 • 10“^ and the confidence bounds 
of the benchmark Monte-Carlo setting do not exceed 6.783 • 10“^. 

For a run-time comparison, we show for different values of M the run-times 
necessary to compute the prices for parameter tupels. Again, the run-times are 
measured for M = 1 and extrapolated for other values of M. Table 15.111 presents 
the results. In Figure [TS] for each M = 1,..., 100 the run-times of the Chebyshev 
interpolation method, including the offline phase, compared to the Monte-Carlo 
method are presented. We observe that for M = 15 both lines intersect and for 
M > 15 the Chebyshev method outperforms its benchmark. Contrary to the case of 
varying only one parameter, the intersection of both lines occurs at a significantly 
lower value of M due to the fact that for each M pricing for parameter tupels 
is performed. 

Additionally, Table 15.111 highlights that for a total number of 50^ parameter 
tuples, the Chebyshev method exhibits a significant decrease in (total) pricing run¬ 
times. For the maximal number of 100^ parameter tuples that we investigated, 
pricing in either model resulted in more than 97% of run-time savings in our im¬ 
plementation. While the computation of 100^ Heston prices via the Monte-Carlo 
method consumes up to 39 days, the Chebyshev method computes the very same 
prices in 23 hours only. Note that only 7 seconds of this time span are consumed by 
actual pricing during the online phase. 
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Heston 


M 

1 

10 

50 

100 

^Cheby / x 
^online ) 


7.1-10-2 

1.8 

7.1 

T^Cheby / x 

offline+online 

8.2-10^ 

8.2-10"^ 

8.2-10^ 

8.2-10"^ 

/j^Monte-Carlo 

3.4-102 

3.4-10^ 

8.4-10^ 

3.4-10® 

^-tCheby 
^ ofiiine+online 
2^Monte-Carlo 

24313.9% 

243.1% 

9.7% 

2.4% 


Table 5.11: Efficiency study for a multivariate lookback option in the Heston model based on 5 
underlyings. Here, we vary two model parameters and compare the Chebyshev results to Monte- 
Carlo. Both methods have been set up to deliver comparable accuracies. With increasing number 
of computed prices, the Chebyshev algorithm increasingly prohts from the initial investment of the 
offline phase. 



M (evaluation of different parameter pairs) 


Figure 5.8: Effiency study for a multivariate lookback option in the Heston model based on 5 
underlyings, varying the two model parameters a and p. Comparison of run-times for Monte- 
Carlo pricing and Chebyshev pricing including the offline phase. Both methods have been set up to 
deliver comparable accuracies. We observe that the Monte-Carlo and the Chebyshev curves roughly 
intersect at M = 15. 


41 














6 Conclusion and Outlook 


This article focuses on applying the Chebyshev method to European option pricing. 
Within this scope, Sections [2H1] establish theoretical convergence results and numer¬ 
ical case studies in Section [5] confirm these findings. Moreover, in experiments with 
Fourier pricing an accuracy of 10“^ is observed with less than ten Chebyshev nodes 
in each parameter, see Figure 15.31 The financial implications of the high precision 
achieved with such a small number of interpolation nodes are twofold. First, it shows 
that there are interesting cases in which we observe an accuracy in the range of ma¬ 
chine precision. In this comfortable situation without methodological risk we can 
ignore the fact that approximations are implemented. Second, compared with other 
sources of risk, already errors from a much lesser accuracy level can be ignored. If 
we agree that an accuracy of 10“^ is satisfactory, already 36-49 interpolation nodes 
for the approximation of call option prices as reported in Figure [5l3] are sufficient. 

Additionally, also the numerical experiments for American, barrier and lookback 
options display promising results. A theoretical error analysis for nonlinear pricing 
problems is beyond the scope of the present article, while we are convinced that fur¬ 
ther investigations in this direction are valuable. For instance our analytic approach 
based on examinations of the Fourier representations can be adopted to barrier op- 
tions in Levy models leading to the involvement of Wiener-Hopf factorizations, see 


Eberlein et al.l (120111 1. In general we expect the regularity analysis to become more 


ch allenging in the pr esence of nonlinearities. For American options the current work 
of lTeichmannl (j2015l l may lead to regularity assertions for American options that are 
inherited from their corresponding European counterparts. 

The theoretical and experimental results of our case studies show that the method 
can perform considerably well when few parameters are varied. As a consequence, 
we recommend the interpolation method for this case and also when solely the 
strike of a plain vanilla option is varied and fast Fourier methods are available. For 
calibration purposes for example, strikes are not given in a discrete logarithmic scale, 
which makes an additional interpolation necessary in order to apply FFT. Here, 
Chebyshev polynomials offer an attractive alternative. In particular, the maturity 
can be used as supplementary free variable. Moreover, for models with a low number 
of parameters, another path could be beneficial: Interpolating the objective function 
of the parameters directly. Then the optimization would boil down to a minimization 
of a tenso rized polynomial, wh ich could be exploited in further research. As may be 


seen from Armenti et ah ( 20151 1. where the present article is applied for the first time. 


this advantage can also be exploited for other optimization procedures in hnance for 
example in risk allocation. 

The multivariate construction of the interpolation and the theoretical error anal¬ 
ysis suggest that the empirically observed error behaviour extends to three and more 
varying parameters, as well, as long as analyticity is provided. More precisely, in 
case of analyticity, the rate is of order p ^ for some constant p depending on 
the domain of analyticity and N the total number of interpolation nodes. For mul¬ 
tivariate polynomial interpolation, the introduction of sparsity techniques promise 
hig her efficiency, as for inst ance by compression techniques for tensors as reviewed 
by Kolda and Bader ( 2009l l. We address the issue of curse of dimensionality fur¬ 
ther in Gafi, Glau and Mair (2015), where we take a different route by replacing the 


42 






















Chebyshev interpolation with an empirical interpolation for Fourier pricing methods. 


A Proof of Theorem 12.2 


Proof. In Sauter and SchwabI ( 2004 . Proof of Lemma 7.3.3) the proof is given for 
the following error bound: 


max 

pS'P 


|/-%(/)|</D2f+V(.;';(i 


^min) 


D_ 

2 


where N is the number of interpolation points in each of the D dimensions, ^min := 
min^i Qi and V the bound of / on B{V, q) with V = [—1,1]^. Here, we extend the 
proof by incorporating the different values of iVj, i = 1,..., L), as well as expressing 
the error bound with the different Qi, i = 1,... , D. 

In general we work with a parameter space V of hyperrectangular structure, 
^ ^ ^ introduced linear transformation in Section 

o we have a transformation r-p : [—1,1]-^ —)• P with 

rv{p) = (pi + ~" 

Let p I—)• Price^ be a function on V. We set Price^ = Price^ o T-p{p). Further, let 
Ijj{Price ){p) be the Chebyshev interpolation of Price^ on [—1,1]^. Then it holds 

Ij^{Price^'^){p) = ij^{Price ^)(-) o rp^(p). 


Therewith, it directly follows 


Price^ — I-^{Price^'^){p) = ^Price — I-^{Price ^)(-)^ o Tj,^{p). 


Applying the error estimation from ISauter and SchwabI (j2004l . Lemma 7.3.3) results 


Price - %(Price('))(-)|^o(p) = |Pnce - /^(Pnce('))(-)|^o([_i,i]i3) 


where V = suppgp(j_]^ Price^, V = sup^g^j-p^^) Price^. Summarizing, the 
transformation rp : [—1,1]^ —?> V does not affect the error analysis, only by ap¬ 
plying the transformation as described in Section 12.11 


(A.l) B{V, q) := B{[p^,p^], Qi) X ...X B{[p^,p^], qd), 

with B{\p,p], g) := r[p^p] o H([—1,1], ^p). Note that Qi is not the radius of the ellipse 
B{\p.,p^], Qi) but of the normed ellipse B{[—1, 1], Qi). Therefore, in the following it 
suffices to show the proof for V = [—1,1]^. 
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As in ISauter and SchwabI (|2nn4l . Proof of Lemma 7.3.3) we introduce the scalar 

iHiipt 


product 


and the Hilbert 




/f \ f 

{f,9)s ■= / —-7- 

nf=i vii - 


:dz 




space 


\J kJ^€A<V-'V-' 

L^{B{V,q)) ■■={/: / is analytic in B{V,q) and ||/||^ := {f,f)g < oo}. 

Following Sauter and Schw^ ( 20041 . Proof of Lemma 7.3.3), we define a complete 
orthonormal system for L?‘{B{V, g)) w.r.t. the scalar product {■,-)q by the scaled 
Chebyshev polynomials 


— D 

f^{z) := c^T^{z) with := , for all jj. G N^. 

Then, for any arbitrary bounded functional E on L^{B{V, g)) we have 
(A.2) |i?(/)|<||i^||,||/||„ 


(A.2) |i?(/)|<||i^||,||/||„ 

where ||F^||g denotes the operator norm. Due to the orthonormality of 
follows that 


/.SN^ 


it 


l|E|l,= sup = 


E i®(b 


^ ' ' ' ... - V 

In the following let E be the error of the Chebyshev polynomial interpolation at a 
fix p G P, 

E{f) :=|/(p)-%(/(-))(p)|. 

Starting with (IA.2I1 . we first focus on ||P||g, 

\\E\\i = Y, = E 

AtSNf /iSNf 

At this step we apply Lemma lA.ll and obtain 


<.u - E 4i^K.)p< E 


E 




4c2 


Overall, using (nf=i 8 ^' + ^) - (nf=i 8 ^/') = U.f=i 8j for 

No and j = 1,..., D and this leads to 


X > 0, /ij G 


I|£|IS<4 E 

, 3 i:fj,i>Ni 


c' <4 


D D 


r, 


< 4 


TT 


D D 


E. 


,-2Ni 




! E E n. 

*=i 



















From this point on we use the convergence of the geometric series since \Qj‘^\ < 

1, j = 1,... ,F), 


D D 


CO oo D 


1^11?S4(i) 

,Ati=0 fijj=0j=l 


i=l 
D D 


CO D—1 CO 

—2fij \ ^ 


Hi) Eft“ME- E 

,/ii=0 Md-i= 0 j=l /^D=0 


Z=1 


4 0)'’E.r- 

^ ^ Z=1 


CO D—1 


= ... = 4 


E E n 

^Ati=0 I1D-1=0 j=l 

D 


2iij 1 


-2 


1 Qd 


D D D 

E^r^niw 


-2 ' 


i=l j=l 

Recalling (IA.2I1 . we have to estimate WfWg, 


L 


mfiz) 


nf=i a/ii- 


:d0< sup |/(Z)| ||1| 

\z&B(v,e) J 


From TT 2 To = 1 it directly follows that ||l||^=(^ 7 r 2 j ||ro||g = vr^ and therewith 

l<7r^- 


Combining the results leads to 


|£(/)l = l/(p) - %(/(-))(p)| < ■ r" ■ 4 (° E n 

i=l ^ 


-2 


D 


Z=1 

D 




-2 


i=l j=l 


□ 

The following lemma shows that the Chebyshev interpolation of a polynomial 
with a degree as most as high as the degree of the interpolating Chebyshev polyno¬ 
mial is exact and furthermore determines an upper bound for interpolating Cheby¬ 
shev polynomials with a higher degree. 

Lemma A.l. For x E [—1,1]'^ it holds 

(A.3) \Tf,{x) - = 0 : Hi<Ni,i = l,...,D, 

(A.4) |r^(x)-%(r^(-))(x)| <2 :3ie{l,...,D}:^ii>Ni. 


45 











Proof. Uniqueness properties o f the Chebyshev int e rpola tion directly imply ()A.3h . 
The proof of (|A.4D is similar to Sauter and Schwab (2004, Proof of Hilfssatz 7.3.1). 
They use the zeros of the Chebyshev polynomial as interpolation points, whereas we 
use the extreme points and therefore, we use a different orthogonality property in this 
proof. We first focus on the one-dimensional case. Recalling (12.2p . the Chebyshev 
interpolation of T^, fj. > N, is given as 


N 


In{T^){x) = '^CjTj{x) with Cj 
j=o 


2^0<j<N 


N 


N 


'^''T^{xk)Tj{xk), j<N, 


k=0 


where Xk den otes th e fc-th extremum of Tat. Here, we can apply the following 
orthogonality (Rivlin (199C, p.54)), 

(A.5) 


N 


'^"T^{xk)Tj{xk) 


= 


0, IJL + j mod (2A^) and |/i — j\ / 0 mod {2N), 

N, + j = 0 mod (2A^) and |/r — j\ =0 mod {2N), 

Y, ^ + J = 0 mod {2N) and |/r — j\ / 0 mod {2N), 

+ J / 0 {2N) and |/r — j\ =0 mod {2N). 

For j < N and fj, > N this yields the existence of 7 < A^ such that 


A:=0 


(A. 6 ) iNiTf,) = T^. 

(ira follows elementarily from the case that for any fj, > N only for one 0 < j < A^ 
the orthogonality can lead to a coefficient cj > 0 . 

Proving the claim, we distinguish several cases. In all of these cases we assume 
that there exists 0 < j < A" such that Ylik=o"'^r-i^k)Tj{xk) 7 ^ 0. We will then show 
that for all other 0 < z < A,z 7 ^ j it follows Ylk=o'''^r-i^k)Tj{xk) = 0. 

Firstly, assume there exists j such that fi + j = 0 mod (2A) and fJ- — j = 
0 mod (2A). Then it directly follows for all 0 < z < A, z 7 ^ j that fj, + i 
0 mod (2A) and /U — z 7 ^ 0 mod (2A). 

Secondly, assume there exists j such that fi + j = 0 mod (2A) and fJ- — j 
0 mod (2A). Analogously, for all 0 < z < A, z 7 ^ j we have fi + i ^ 0 mod (2A) and 
additionally from fj. + j = 0 mod (2 A) it follows that fJ. + j — 2 A = 0 mod (2 A) and 
therewith for all 0 < z < A, z 7 ^ j we have ^ — i>^ + j — 2N which is equivalent 
to /z — z 7 ^ 0 mod (2A). 

Similar argumentation holds for the third case n+j 7 ^ 0 mod (2A) and |/z —= 
0 mod (2A). 

Therewith, (IA. 6 D holds and it directly follows that |r^ — lN{Tfj,)\ < \T^\ + 
|/ 7 v(T^)| < 1 -|- 1 = 2. Thus (IA.4I1 holds in the one-dimensional case. The extension 
to the multivariate case follows analogously by applying the triangle inequality and 
inserting the one-dimensional result to each tensor component. □ 
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